{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-point stress approximation (MPSA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PorePy includes a multi-point stress approximation discretization for the linear elasticity problem:\n",
    "\\begin{equation}\n",
    "\\nabla\\cdot \\sigma = -\\vec f,\\quad \\vec x \\in \\Omega\n",
    "\\end{equation}\n",
    "where $\\vec f$ is a body force, and the stress $\\sigma$ is given as a linear function of the displacement\n",
    "\\begin{equation}\n",
    "\\sigma = C: \\nabla \\vec u.\n",
    "\\end{equation}\n",
    "\n",
    "The convention in PorePy is that tension is positive. This means that the Cartesian component of the traction $\\vec T = \\sigma \\cdot \\vec n$, for a direction $\\vec r$ is positive number if the inner product $\\vec T\\cdot \\vec r$ is positive. The displacements will give the difference between the initial state of the rock and the deformed state. If we consider a point in its initial state $\\vec x \\in \\Omega$ and let $\\vec x^* \\in \\Omega$ be the same point in the deformed state, to be consistent with the convention we used for traction, the displacements are given by $\\vec u = \\vec x^* - \\vec x$, that is, $u$ points from the initial state to the finial state.\n",
    "\n",
    "To close the system we also need to define a set of boundary conditions. Here we have three posibilities, Neumman conditions, Dirichlet conditions or Robin conditions, and we divide the boundary into three disjont sets $\\Gamma_N$, $\\Gamma_D$ and $\\Gamma_R$ for the three different types of boundary conditions\n",
    "\\begin{equation}\n",
    "\\vec u = g_D, \\quad \\vec x \\in \\Gamma_D\\\\\n",
    "\\sigma\\cdot n = g_N, \\quad \\vec x \\in \\Gamma_N\\\\\n",
    "\\sigma\\cdot n + W \\vec u= g_R,\\quad \\vec x\\in \\Gamma_R.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve this system we first have to create the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import porepy as pp\n",
    "\n",
    "# Create grid\n",
    "n = 5\n",
    "g = pp.CartGrid([n, n])\n",
    "g.compute_geometry()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to define the stress tensor $C$. In PorePy, the constitutive law\n",
    "\\begin{equation}\n",
    "\\sigma = C:\\nabla u = 2  \\mu  \\epsilon +\\lambda  \\text{trace}(\\epsilon) I, \\quad \\epsilon = \\frac{1}{2}(\\nabla u + (\\nabla u)^\\top)\n",
    "\\end{equation}\n",
    "is implemented, and to get the tensor for this law we call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create stiffness matrix\n",
    "lam = np.ones(g.num_cells)\n",
    "mu = np.ones(g.num_cells)\n",
    "C = pp.FourthOrderTensor(mu, lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we need to define boundary conditions. We set the bottom boundary as a Dirichlet boundary, and the other boundaries are set to Neuman."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define boundary type\n",
    "dirich = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))\n",
    "bound = pp.BoundaryConditionVectorial(g, dirich, [\"dir\"] * dirich.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize the stresses by using the multi-point stress approximation (for details, please see: E. Keilegavlen and J. M. Nordbotten. “Finite volume methods for elasticity with weak symmetry”. In: International Journal for Numerical Methods in Engineering (2017))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the values we put on the boundaries. We clamp the bottom boundary, and push down by a unitary traction on the top boundary. Note that the value of the Neumann condition given on a face $\\pi$ is the integrated traction $\\int_\\pi g_N d\\vec x$, hence the multiplication by face areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_faces = np.ravel(np.argwhere(g.face_centers[1] > n - 1e-10))\n",
    "bot_faces = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))\n",
    "\n",
    "u_b = np.zeros((g.dim, g.num_faces))\n",
    "u_b[1, top_faces] = -1 * g.face_areas[top_faces]\n",
    "u_b[:, bot_faces] = 0\n",
    "\n",
    "u_b = u_b.ravel(\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize this system using the `Mpsa` class. We assume zero body forces $f=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameter_keyword = \"mechanics\"\n",
    "\n",
    "mpsa_class = pp.Mpsa(parameter_keyword)\n",
    "f = np.zeros(g.dim * g.num_cells)\n",
    "\n",
    "specified_parameters = {\n",
    "    \"fourth_order_tensor\": C,\n",
    "    \"source\": f,\n",
    "    \"bc\": bound,\n",
    "    \"bc_values\": u_b,\n",
    "}\n",
    "data = pp.initialize_default_data(g, {}, parameter_keyword, specified_parameters)\n",
    "mpsa_class.discretize(g, data)\n",
    "A, b = mpsa_class.assemble_matrix_rhs(g, data)\n",
    "\n",
    "u = np.linalg.solve(A.toarray(), b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unknowns are ordered so that u[0] and u[1] contain the displacement in the x- and y-direction in cell 0, respectively, u[2] gives the x-displacement in cell 1 etc. Thus we can plot the y component of the displacement by writing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh0AAAHHCAYAAAAbLeozAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3hUlEQVR4nO3df1iUdb7/8Rc/AhQdWBRFNhTNEq3U1ORoVpZUSuuqp2PlYUtdV1tXKtOt6LubaeWSnVrth4t1ThttV1z9sNVtrbVDmlqmgri01kFSDwWLIhbBBCbCzHz/MOZEDMpw38w9Mz4f13VfNPfcP95z77Tz6vPjvkNcLpdLAAAAXSzU6gIAAMC5gdABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABoMt8/vnnCgkJUW5ubqf2z83NVUhIiD7//HNT6wJgDUIH8D2HDx/WHXfcoUGDBikqKko2m01XXHGFnnrqKX377bdWl+d25MgRLV++XMXFxV7tV1ZWpszMTF100UXq3r27unfvrmHDhmnRokX6xz/+0TXFAsB3wq0uAPAXb7/9tmbOnKnIyEjdfvvtuuSSS3Tq1Cl9+OGHuvfee/Xpp5/q+eeft7pMSadDx4oVK5ScnKyRI0d2aJ9NmzbplltuUXh4uDIyMjRixAiFhobqwIED+vOf/6ycnByVlZVpwIABptU5YMAAffvttzrvvPNMOyaAwEXoAHS6BeDWW2/VgAEDtHXrVvXr18/93qJFi3To0CG9/fbbPq2poaFB0dHRphzr8OHD7s+3ZcuWVp9PklatWqU//OEPCg31vvHzTHWGhIQoKiqqUzUDCD50rwCSHn/8cdXX1+uFF15o84MsSYMHD9bdd999xmN88MEHmjlzpvr376/IyEglJSXpnnvu6VC3TMvYhe3bt+tXv/qV+vTpo/PPP9/jttu2bdPll18uSZo7d65CQkLOOm7i8ccfV0NDg1588UWPny88PFx33XWXkpKSTKtT8m5Mx6effqprr71W3bp10/nnn69HH31UTqfzrPsBCBy0dACS/vrXv2rQoEEaP358p4/xxhtv6MSJE1q4cKF69eqlgoICPfPMM/rnP/+pN954o0PH+NWvfqX4+HgtW7ZMDQ0NHrcZOnSoHn74YS1btkwLFizQlVdeKUlnrH3Tpk0aPHiwUlNTvf9gnazTG1VVVbrmmmvU3NysrKwsRUdH6/nnn1e3bt1MqBaAvyB04Jxnt9tVWVmpadOmGTrOqlWrWv1ILliwQIMHD9b/+3//T+Xl5erfv/9ZjxEXF6ctW7YoLCys3W369u2rKVOmaNmyZRo3bpx+9rOfnfGYdrtdR44c0fTp09u8V1tbq+bmZvfr6OjoDv3Qd6ROb6xatUrHjx/Xnj17NHbsWEnS7NmzdeGFF5pyfAD+ge4VnPPsdrskqWfPnoaO8/0f64aGBn355ZcaP368XC6X/v73v3foGPPnzzfth7xFy+fr0aNHm/cmTpyo+Ph497J27VpL6nznnXf0L//yL+7AIUnx8fHKyMgw7RwArEdLB855NptNkvTNN98YOk55ebmWLVumt956S19//XWr9+rq6jp0jIEDBxqqwZOWMFVfX9/mveeee07ffPONjh07dtYWk+8zu84vvvjCY9fPkCFDTD0PAGsROnDOs9lsSkxM1CeffNLpYzgcDl133XWqqanR/fffr5SUFEVHR6uyslJz5szp8IDIrhjDEBMTo379+nn8fC0/9N7efIuxFgA6g+4VQNJPfvITHT58WLt27erU/vv379dnn32mJ598Uvfff7+mTZumtLQ0JSYmmlzpaSEhIV5tf+ONN+rQoUMqKCjoknqMGjBggA4ePNhmfWlpqQXVAOgqhA5A0n333afo6Gj94he/0LFjx9q8f/jwYT311FPt7t8yvsHlcrnXuVyuM+5jRMt9MWprazu0/X333afu3bvr5z//ucfP9/26rZCenq7du3e3CkXHjx/XK6+8YmFVAMxG9wog6YILLlBeXp5uueUWDR06tNUdST/66CO98cYbmjNnTrv7p6Sk6IILLtCvf/1rVVZWymaz6c0332wztsPMemNjY7Vu3Tr17NlT0dHRSk1NbXesxYUXXqi8vDzNmjVLQ4YMcd+R1OVyqaysTHl5eQoNDT3jPTe60n333aeXX35ZkydP1t133+2eMjtgwABuzw4EExcAt88++8w1f/58V3JysisiIsLVs2dP1xVXXOF65plnXCdPnjzjvv/zP//jSktLc/Xo0cPVu3dv1/z5810ff/yxS5LrxRdfPOO+L774okuSq7CwsMO1/uUvf3ENGzbMFR4e3qFzuFwu16FDh1wLFy50DR482BUVFeXq1q2bKyUlxfXLX/7SVVxcfNb9va2zrKysw7X94x//cF199dWuqKgo149//GPXI4884nrhhRdcklxlZWUdOh8A/xbiclncrgoAAM4JjOkAAAA+QegAAAA+QegAAAA+YWnoWL58ufsJmS1LSkqKlSUBAIAuYvmU2Ysvvljvvfee+3V4uOUlAQCALmD5L3x4eLgSEhKsLgMAAHQxy0PHwYMHlZiYqKioKI0bN07Z2dntPgK8sbFRjY2N7tdOp1M1NTXq1auX17eFBgCcW1wul7755hslJiYqNLTrRhecPHlSp06dMnyciIgIRUVFmVCR/7D0Ph1/+9vfVF9fryFDhujo0aNasWKFKisr9cknn3h8zPjy5cu1YsUKCyoFAASLioqKLrv77smTJxXfrZvaPtPZewkJCSorKwuq4OFXNwerra3VgAED9Pvf/17z5s1r8/4PWzrq6uq+axWZLIkuGu8dlvSBpGmS4i2uJdAclLRd0k2SeltcS6A5KGmr+N51Vst3b6qkXhbXEmiqJG1WbW2tYmJiuuQMdrtdMTExukdSpIHjNEpardO/czabzZzi/IDl3SvfFxsbq4suukiHDh3y+H5kZKQiIz39z5ggKbkrSwtS9u/+/lhS1zwNNXjVfff3x98t6LjvXzu+d95ruX79xPXrHF90x0dLMtI+4Vc/zibyq/t01NfX6/Dhw+rXr5/VpQAA0GnnmbAEI0tDx69//Wtt375dn3/+uT766CPNmDFDYWFhmjVrlpVlAQBgSLgJS1epqalRRkaGbDabYmNjNW/ePNXXtz8KpaamRnfeeaeGDBmibt26qX///rrrrrtUV1fX7j7tsbQF55///KdmzZqlr776SvHx8ZowYYJ2796t+Hj6eQEA6AoZGRk6evSo8vPz1dTUpLlz52rBggXKy8vzuP2RI0d05MgRPfHEExo2bJi++OIL/fKXv9SRI0e0fv16r85taeh49dVXrTw9AABdIlzGukiazSrkB0pKSrR582YVFhZqzJgxkqRnnnlG6enpeuKJJ5SY2Hac0CWXXKI333zT/fqCCy7QypUr9bOf/UzNzc1e3dTTr8Z0AAAQDPy1e2XXrl2KjY11Bw5JSktLU2hoqPbs2dPh47TMqvH2LuLBOkAWAICAZ7fbW71ufxZnx1RVValPnz6t1oWHhysuLk5VVVUdOsaXX36pRx55RAsWLPD6/LR0AABgMrNmryQlJSkmJsa9ZGdnezxfVlZWmweo/nA5cOCA4c9lt9t14403atiwYVq+fLnX+9PSAQCAyYx2kbTsW1FR0ermYO21cixdulRz5sw54zEHDRqkhIQEVVdXt1rf3Nysmpqasz4H7ZtvvtHkyZPVs2dPbdiwQeed5/2oFUIHAAB+ymazdeiOpPHx8R2a+Tlu3DjV1taqqKhIo0ePliRt3bpVTqdTqamp7e5nt9t1ww03KDIyUm+99Vanb81O9woAACZrmb3S2aWrWgSGDh2qyZMna/78+SooKNDOnTuVmZmpW2+91T1zpbKyUikpKSooKJB0OnBcf/31amho0AsvvCC73a6qqipVVVXJ4XB4dX5aOgAAMJlZ3Std4ZVXXlFmZqYmTZqk0NBQ3XTTTXr66afd7zc1Nam0tFQnTpyQJO3bt889s2Xw4MGtjlVWVqbk5OQOn5vQAQDAOSQuLq7dG4FJUnJysr7/LNiJEyfKrGfDEjoAADCZ0eenBOuzVwgdAACYjNDhGaEDAACT+fOYDisxewUAAPhEsIYpAAAsY/SBb8H64xysnwsAAMvQveIZ3SsAAMAngjVMAQBgGWaveEboAADAZHSveEb3CgAA8IlgDVMAAFiG2SueBevnAgDAMnSveEb3CgAA8IlgDVMAAFiG2SueEToAADAZ3SueBevnAgDAMgwk9YwxHQAAwCeCNUwBAGAZxnR4RugAAMBkjOnwjO4VAADgE8EapgAAsEx4mHReiIH9XZIcppXjNwgdAACYLDxcCid0tEH3CgAA8AlaOgAAMNl5BrtXznOZV4s/IXQAAGAyU7pXghDdKwAAwCdo6QAAwGTnhUnnGfjP+vOc5tXiTwgdAACYLUzG+hIMdM34M0IHAABmC5ex0BGkLR2M6QAAAD5BSwcAAGajpcMjQgcAAGYjdHhE9woAAPAJWjoAADBbqE7PYEErhA4AAMwWLmOhI0inzNK9AgAAfIKWDgAAzEZLh0eEDgAAzBYmxnR4QPcKAADnkJqaGmVkZMhmsyk2Nlbz5s1TfX19h/Z1uVyaMmWKQkJCtHHjRq/PTegAAMBs4SYsXSQjI0Offvqp8vPztWnTJu3YsUMLFizo0L5r1qxRSEjn+37oXgEAwGxh8stf2JKSEm3evFmFhYUaM2aMJOmZZ55Renq6nnjiCSUmJra7b3FxsZ588knt3btX/fr169T5aekAAMBsYSYsXWDXrl2KjY11Bw5JSktLU2hoqPbs2dPufidOnNC///u/a+3atUpISOj0+f0whwEAAEmy2+2tXkdGRioyMrLTx6uqqlKfPn1arQsPD1dcXJyqqqra3e+ee+7R+PHjNW3atE6fW6KlAwAA85k0piMpKUkxMTHuJTs72+PpsrKyFBIScsblwIEDnfoob731lrZu3ao1a9Z0av/vo6UDAACzmTQYtKKiQjabzf26vVaOpUuXas6cOWc81qBBg5SQkKDq6upW65ubm1VTU9Nut8nWrVt1+PBhxcbGtlp/00036corr9S2bdvO+jlaEDoAAPBTNputVehoT3x8vOLj48+63bhx41RbW6uioiKNHj1a0ulQ4XQ6lZqa6nGfrKws/eIXv2i17tJLL9Xq1as1derUDnyK/0PoAADAbF087bWzhg4dqsmTJ2v+/Plat26dmpqalJmZqVtvvdU9c6WyslKTJk3Sn/70J40dO1YJCQkeW0H69++vgQMHenV+xnQAAGC2lqfMdnbpwl/nV155RSkpKZo0aZLS09M1YcIEPf/88+73m5qaVFpaqhMnTph+bj/MYQAAoKvExcUpLy+v3feTk5PlcrnOeIyzvd8eQgcAAGYz2r3Sud90v0foAADAbIQOjxjTAQAAfIKWDgAAzGb0VuZOswrxL4QOAADMRveKR4QOAADMZvQps0Ha0uE3Yzoee+wxhYSEaPHixVaXAgAAuoBftHQUFhbqueee0/Dhw60uBQAA44yO6eiiR9tbzfKWjvr6emVkZOg///M/9aMf/cjqcgAAMM6kp8wGG8s/1qJFi3TjjTcqLS1Njz766Bm3bWxsVGNjo/u13W7/7p++khTRdUUGra+/+1t9xq3gSc13f7l23uPaGdNy/b60tIrA9JXVBZzzLA0dr776qvbt26fCwsIObZ+dna0VK1Z4eOev5hZ2zllvdQEB7FWrCwhgfO+M+bPVBeBMjLZWBOlAUstCR0VFhe6++27l5+crKiqqQ/s88MADWrJkifu13W5XUlKSNOE3Uu+hXVVq8KrYKRXlSGmPSD/y7kmB57wvdkoFOdJ1j0hxXDuvfL5T2sP3rtNavnsTH5FiuX5e+bJE+nClb85F6PDIstBRVFSk6upqjRo1yr3O4XBox44devbZZ9XY2KiwsNYjaSIjIxUZGdn2YIOulwZc1dUlB6eiHOmidOnHo86+LVoryJFSuHadsofvnSEFOdLgdKkf188rX+zwXeiAR5aFjkmTJmn//v2t1s2dO1cpKSm6//772wQOAAACRsuj7Y3sH4QsCx09e/bUJZdc0mpddHS0evXq1WY9AAABxWj3isOsQvxLkGYpAADgbyyfMvt927Zts7oEAACMo6XDI78KHQAABAXuSOoRoQMAALPR0uERYzoAAIBP0NIBAIDZjD7avtmsQvwLoQMAALMZ7V4J0l9nulcAAIBPBGmWAgDAQsxe8YjQAQCA2ehe8YjuFQAA4BNBmqUAALAQLR0eBenHAgDAQjxl1qMg/VgAAMDf0NIBAIDZ6F7xKEg/FgAAFiJ0eBSkHwsAAAtxnw6PGNMBAAB8gpYOAADMRveKR0H6sQAAsJDRp8zSvQIAANB5tHQAAGA2ulc8CtKPBQCAhZi94hHdKwAAnENqamqUkZEhm82m2NhYzZs3T/X19Wfdb9euXbr22msVHR0tm82mq666St9++61X5yZ0AABgtnATli6SkZGhTz/9VPn5+dq0aZN27NihBQsWnHGfXbt2afLkybr++utVUFCgwsJCZWZmKjTUuxhB9woAAGbz0zEdJSUl2rx5swoLCzVmzBhJ0jPPPKP09HQ98cQTSkxM9LjfPffco7vuuktZWVnudUOGDPH6/LR0AABwjti1a5diY2PdgUOS0tLSFBoaqj179njcp7q6Wnv27FGfPn00fvx49e3bV1dffbU+/PBDr89P6AAAwGwtj7bv7PLdr7Pdbm+1NDY2GiqrqqpKffr0abUuPDxccXFxqqqq8rjP//7v/0qSli9frvnz52vz5s0aNWqUJk2apIMHD3p1fkIHAABmM2lMR1JSkmJiYtxLdna2x9NlZWUpJCTkjMuBAwc69VGcTqck6Y477tDcuXN12WWXafXq1RoyZIj++Mc/enUsxnQAAGA2k8Z0VFRUyGazuVdHRkZ63Hzp0qWaM2fOGQ85aNAgJSQkqLq6utX65uZm1dTUKCEhweN+/fr1kyQNGzas1fqhQ4eqvLz8jOf8IUIHAAB+ymaztQod7YmPj1d8fPxZtxs3bpxqa2tVVFSk0aNHS5K2bt0qp9Op1NRUj/skJycrMTFRpaWlrdZ/9tlnmjJlSgc+xf+hewUAALMZGc9h9MZiZzB06FBNnjxZ8+fPV0FBgXbu3KnMzEzdeuut7pkrlZWVSklJUUFBgSQpJCRE9957r55++mmtX79ehw4d0oMPPqgDBw5o3rx5Xp2flg4AAMzmp1NmJemVV15RZmamJk2apNDQUN100016+umn3e83NTWptLRUJ06ccK9bvHixTp48qXvuuUc1NTUaMWKE8vPzdcEFF3h1bkIHAADnkLi4OOXl5bX7fnJyslwuV5v1WVlZre7T0RmEDgAAzMaj7T0idAAAYDY/7l6xEgNJAQCATwRplgIAwEI82t4jQgcAAGaje8UjulcAAIBPBGmWAgDAQsxe8YjQAQCA2RjT4RGhAwAAszGmwyPGdAAAAJ8I0iwFAICFaOnwKEg/FgAAFiJ0eET3CgAA8IkgzVIAAFjHFSq5DMxAcQVpkwChAwAAkznCTy9G9g9GQZqlAACAvwnSLAUAgHVo6fAsSD8WAADWaQ4LUXNYiIH9XZJc5hXkJ+heAQAAPkFLBwAAJnOEh8sR3vmWDke4S1KTeQX5CUIHAAAmc4SFyWGge8URRugAAAAd4FSYHOp86HAG4XgOiTEdAADAR2jpAADAZM0KU7OBlo7mIG3pIHQAAGAyh8LkMNCZ4JDTxGr8B90rAADAJ2jpAADAZMZbOjrfNePPCB0AAJiM0OEZ3SsAAMAnaOkAAMBktHR4RugAAMBkDoWpmdDRhqXdKzk5ORo+fLhsNptsNpvGjRunv/3tb1aWBACAYQ6FG16CkaWh4/zzz9djjz2moqIi7d27V9dee62mTZumTz/91MqyAABAF7A0Sk2dOrXV65UrVyonJ0e7d+/WxRdfbFFVAAAY41CoHAozsH9w8pv2G4fDoTfeeEMNDQ0aN26cx20aGxvV2Njofm2320//w1elUkQPX5QZXGrLTv89XmJtHYHo6++uXTXXzms1fO8Mafnufcn189pXpT471emBpISOHwpxuVyW3uB9//79GjdunE6ePKkePXooLy9P6enpHrddvny5VqxY4eMKAQDBpK6uTjabrUuObbfbFRMTox11F6mHrfOho97u0FUxn3VprVawPHScOnVK5eXlqqur0/r16/Vf//Vf2r59u4YNG9ZmW08tHUlJSUr9zdWKGxrvy7KDwpGdX+jjnEJd88iV+tHAGKvLCSjlOyu1N+fvmvTIFfrRwOD5PwRf+GLnERXmfMz3rpNavntXPDJJtoE/srqcgFJTclx7Vm73SejYWjfUcOi4NqYk6EKH5d0rERERGjx4sCRp9OjRKiws1FNPPaXnnnuuzbaRkZGKjIxss37A9YOVdFVyV5calD7OKdSF6YOUOCrB6lICzt6cv+ui9IFKHNXX6lICTmHOx3zvDNib83cNTL9IfUclWl1KQKnY8bn2rNzuk3M5FW6oe8XJlFnfcDqdrVozAACAeWpqapSRkSGbzabY2FjNmzdP9fX1Z9ynqqpKt912mxISEhQdHa1Ro0bpzTff9PrclrZ0PPDAA5oyZYr69++vb775Rnl5edq2bZveffddK8sCAMAQfx5ImpGRoaNHjyo/P19NTU2aO3euFixYoLy8vHb3uf3221VbW6u33npLvXv3Vl5enm6++Wbt3btXl112WYfPbWnoqK6u1u23366jR48qJiZGw4cP17vvvqvrrrvOyrIAADDEX0NHSUmJNm/erMLCQo0ZM0aS9Mwzzyg9PV1PPPGEEhM9d9l99NFHysnJ0dixYyVJv/3tb7V69WoVFRUFTuh44YUXrDw9AAB+zX1riO+0N7axo3bt2qXY2Fh34JCktLQ0hYaGas+ePZoxY4bH/caPH6/XXntNN954o2JjY/X666/r5MmTmjhxolfn97sxHQAABLqWm4N1fjn985yUlKSYmBj3kp2dbaiuqqoq9enTp9W68PBwxcXFqaqqqt39Xn/9dTU1NalXr16KjIzUHXfcoQ0bNrgngnSU5bNXAAAINs0KU7OB7pVmnb6bRUVFRasps+21cmRlZWnVqlVnPGZJSedvKPfggw+qtrZW7733nnr37q2NGzfq5ptv1gcffKBLL720w8chdAAAYDKjD21rGdPR8kDUs1m6dKnmzJlzxm0GDRqkhIQEVVdXt1rf3NysmpoaJSR4nsJ++PBhPfvss/rkk0/cjygZMWKEPvjgA61du1br1q07a30tCB0AAAS4+Ph4xcef/SaZ48aNU21trYqKijR69GhJ0tatW+V0OpWamupxnxMnTkiSQkNbj8gICwuT0+n0qk7GdAAAYDKnofEcYXIa6Jo5k6FDh2ry5MmaP3++CgoKtHPnTmVmZurWW291z1yprKxUSkqKCgoKJEkpKSkaPHiw7rjjDhUUFOjw4cN68sknlZ+fr+nTp3t1fkIHAAAmMzaI1Nh027N55ZVXlJKSokmTJik9PV0TJkzQ888/736/qalJpaWl7haO8847T++8847i4+M1depUDR8+XH/605/00ksvtfustPbQvQIAwDkkLi7ujDcCS05O1g8fy3bhhRd26g6kP0ToAADAZM0KNTh7xbuxEoGC0AEAgMmMz16x9AHwXYYxHQAAwCdo6QAAwGTGn71C9woAAOgAQodndK8AAACfoKUDAACTOQw+eyVYWzoIHQAAmIzZK54ROgAAMFnLo+07v7/j7BsFIMZ0AAAAn6ClAwAAkxmfvdJ1z16xEqEDAACTETo8o3sFAAD4BC0dAACYzPiU2eBs6SB0AABgMuNTZoPzPh10rwAAAJ+gpQMAAJMxkNQzQgcAACYzfnOw4OyICM5PBQAA/A4tHQAAmKzZ4OwVI/v6M69bOmbPnq0dO3Z0RS0AAASFltkrRpZg5HXoqKurU1pami688EL97ne/U2VlZVfUBQBAwHJ+N5C0s4uTlo7TNm7cqMrKSi1cuFCvvfaakpOTNWXKFK1fv15NTU1dUSMAAAgCnRpIGh8fryVLlujjjz/Wnj17NHjwYN12221KTEzUPffco4MHD5pdJwAAAcNIK4fR6bb+zNDslaNHjyo/P1/5+fkKCwtTenq69u/fr2HDhmn16tVm1QgAQEBpmTLb+SU4J5d6/amampr05ptv6ic/+YkGDBigN954Q4sXL9aRI0f00ksv6b333tPrr7+uhx9+uCvqBQAAAcrr4bH9+vWT0+nUrFmzVFBQoJEjR7bZ5pprrlFsbKwJ5QEAEHiaFaYwpsy24XXoWL16tWbOnKmoqKh2t4mNjVVZWZmhwgAACFTGH/gWnFNmvf5Ut912W1fUAQAAglxwRikAACzkNDgDJVjv00HoAADAZDxl1rPgnJMDAAD8Di0dAACYrFlhCmX2ShuEDgAATHa6e8XI7BVCBwAA6ADGdHjGmA4AAOATtHQAAGAyWjo8I3QAAGAy7tPhGd0rAACcQ1auXKnx48ere/fuHX5Omsvl0rJly9SvXz9169ZNaWlpOnjwoNfnJnQAAGCyZoUZXrrKqVOnNHPmTC1cuLDD+zz++ON6+umntW7dOu3Zs0fR0dG64YYbdPLkSa/OTfcKAAAmcyhMoX46ZXbFihWSpNzc3A5t73K5tGbNGv32t7/VtGnTJEl/+tOf1LdvX23cuFG33nprh89NSwcAAH7Kbre3WhobG31eQ1lZmaqqqpSWluZeFxMTo9TUVO3atcurYxE6AAAwWcvsFSOLJCUlJSkmJsa9ZGdn+/yzVFVVSZL69u3ban3fvn3d73UU3SsAAJjMYfA26C2ho6KiQjabzb0+MjLS4/ZZWVlatWrVGY9ZUlKilJSUTtdkBkIHAAB+ymaztQod7Vm6dKnmzJlzxm0GDRrUqRoSEhIkSceOHVO/fv3c648dO6aRI0d6dSxCBwAAJmtWmEJ8+MC3+Ph4xcfHd/p8ZzJw4EAlJCRoy5Yt7pBht9u1Z88er2bASIzpAADAdE6Fy2FgcXZhm0B5ebmKi4tVXl4uh8Oh4uJiFRcXq76+3r1NSkqKNmzYIEkKCQnR4sWL9eijj+qtt97S/v37dfvttysxMVHTp0/36ty0dAAAYDKHwZaOrpwyu2zZMr300kvu15dddpkk6f3339fEiRMlSaWlpaqrq3Nvc99996mhoUELFixQbW2tJkyYoM2bNysqKsqrcxM6AAA4h+Tm5p71Hh0ul6vV65CQED388MN6+OGHDZ2b0AEAgMkcCjXY0hGcox8IHQAAmOz0QFDfDSQNFMEZpQAAgN+hpQMAAJM5FK4QQ89eCc6f5+D8VAAAWMj5vVuZd3b/YET3CgAA8AlaOgAAMJnD4EDSrrxPh5UsbenIzs7W5Zdfrp49e6pPnz6aPn26SktLrSwJAADDzHrKbLCxNHRs375dixYt0u7du5Wfn6+mpiZdf/31amhosLIsAADQBSztXtm8eXOr17m5uerTp4+Kiop01VVXWVQVAADGNCtULm4O1oZfjelouc97XFycx/cbGxvV2Njofm232yVJX5d+qYgeEV1fYJCxl30tSfqy5CuLKwk8X5ed/q4eL6mxuJLA83XZ6X9v+d51Tst3r6bkuMWVBJ6vS7/02blOT3llyuwPhbh+eIN1izidTv30pz9VbW2tPvzwQ4/bLF++XCtWrPBxZQCAYFJXVyebzdYlx7bb7YqJidEFdTsVZuvR6eM47PU6HHNFl9ZqBb+JUosWLdInn3zSbuCQpAceeEBLlixxv7bb7UpKStJvZkhDE31RZXDZ+ZmUky89ki4N7GV1NYFl5/9KOTulRyZLAz03zKEdO8uknF187zrL/d27WRoYb3U1gaXkiLRyg9VVnNv8InRkZmZq06ZN2rFjh84///x2t4uMjFRkZGSb9ddfKl01tCsrDF45+VL6MGlUktWVBJ6cnVL6UGlU+19ZtCNnF987I3J2SukjpVEDra4ksOwo8V3ocBqcMhusNwezNHS4XC7deeed2rBhg7Zt26aBA/k3CAAQ+JoVplBCRxuWho5FixYpLy9Pf/nLX9SzZ09VVVVJkmJiYtStWzcrSwMAACazNHTk5ORIkiZOnNhq/Ysvvqg5c+b4viAAAEzgUJhcBn5iaenoAn4ycQYAAFOdDh10r/xQcN59BAAA+B2/mL0CAEAwoaXDM0IHAAAmczjD5HIaCB0G9vVndK8AAACfoKUDAACTOZrD5GzufGuFy8C+/ozQAQCAyRzN4Qpp7vxPrMvAvv4sOD8VAAAWcjSHKsRQS0dwjn4Izk8FAAD8Di0dAACYzNEcZrClgzEdAACgA5qbwxTSROj4IbpXAACAT9DSAQCAyVyOcLkcBn5ijezrx4LzUwEAYKXmsNOLkf2DEN0rAADAJ2jpAADAbLR0eEToAADAbI4QqTnE2P5BiO4VAADgE7R0AABgtubvFiP7ByFCBwAAZiN0eEToAADAbIQOjxjTAQDAOWTlypUaP368unfvrtjY2LNu39TUpPvvv1+XXnqpoqOjlZiYqNtvv11Hjhzx+tyEDgAAzNYsqcnA0oUtHadOndLMmTO1cOHCDm1/4sQJ7du3Tw8++KD27dunP//5zyotLdVPf/pTr89N9woAAGZzfLcY2b+LrFixQpKUm5vboe1jYmKUn5/fat2zzz6rsWPHqry8XP379+/wuQkdAAD4Kbvd3up1ZGSkIiMjLarm/9TV1SkkJKRD3TPfR/cKAABmazZhkZSUlKSYmBj3kp2d7dvP4cHJkyd1//33a9asWbLZbF7tS0sHAABmM2n2SkVFRasf9vZaObKysrRq1aozHrKkpEQpKSkGijo9qPTmm2+Wy+VSTk6O1/sTOgAA8FM2m61DrQlLly7VnDlzzrjNoEGDDNXSEji++OILbd261etWDonQAQCA+Xx8n474+HjFx8cbOOGZtQSOgwcP6v3331evXr06dRzGdAAAYDaHjI3n6MLZK+Xl5SouLlZ5ebkcDoeKi4tVXFys+vp69zYpKSnasGGDpNOB49/+7d+0d+9evfLKK3I4HKqqqlJVVZVOnTrl1blp6QAA4ByybNkyvfTSS+7Xl112mSTp/fff18SJEyVJpaWlqqurkyRVVlbqrbfekiSNHDmy1bG+v09HEDoAADCbH98GPTc396z36HC5XO5/Tk5ObvXaCEIHAABm8+PQYSVCBwAAZmu5nbmR/YMQA0kBAIBP0NIBAIDZ/PjZK1YidAAAYLaWKbNG9g9CdK8AAACfoKUDAACzMXvFI0IHAABmI3R4RPcKAADwCVo6AAAwGy0dHhE6AAAwG7NXPKJ7BQAA+AQtHQAAmI3uFY8IHQAAmK1JUpjB/YMQoQMAALNxG3SPGNMBAAB8gpYOAADMxpgOjwgdAACYjSmzHtG9AgAAfIKWDgAAzNYsY7NX6F4BAAAd0iRjfQlBOmWW7hUAAOATtHQAAGA27tPhEaEDAACzMXvFI7pXAACAT9DSAQCA2Zpl7D/rmb0CAAA6pElSiMH9gxChAwAAszGQ1CPGdAAAAJ+gpQMAALMxpsMjQgcAAGZjyqxHdK8AAACfsDR07NixQ1OnTlViYqJCQkK0ceNGK8sBAMAcTSYsQcjS0NHQ0KARI0Zo7dq1VpYBAIC5HCYsQcjSMR1TpkzRlClTrCwBAAD4SEANJG1sbFRjY6P7td1ulySVHpV6RFlVVeAqO376b8kxa+sIRGVfnf7LtfNeWc3pv1y7znF/9yqtrSMQlR714cmaZezmYMxesV52drZWrFjRZv2C/7SgmCDys5etriBw/SzP6goCF987Y35Gr7R/I3R4FFCh44EHHtCSJUvcr+12u5KSkpQuqZ91ZQWsQ5K2S7pG0o8sriXQlEvaK65dZ3DtjGm5ftMk9ba4lkBzVNI7Vhdxjguo0BEZGanIyMg26y+WdJHvywkK2yVdKCnR6kIC0F5x7TqLa2fMXkmXShpgdSEB5jP5MHQYbanowpaOlStX6u2331ZxcbEiIiJUW1vr1f6//OUv9dxzz2n16tVavHixV/tynw4AAMzmx7NXTp06pZkzZ2rhwoVe77thwwbt3r1biYmd+08GS1s66uvrdejQIffrsrIyFRcXKy4uTv3797ewMgAADPDjlo6WsZG5uble7VdZWak777xT7777rm688cZOndvS0LF3715dc8017tct4zVmz57t9cUAACDYtMzSbNHeMIOu5nQ6ddttt+nee+/VxRdf3OnjWBo6Jk6cKJfLZWUJAACYz6SWjqSkpFarH3roIS1fvtzgwb23atUqhYeH66677jJ0nIAaSAoAQEBolmTkv6m/G9NRUVEhm83mXt1eK0dWVpZWrVp1xkOWlJQoJSXF61KKior01FNPad++fQoJMTIPmNABAIDfstlsrUJHe5YuXao5c+accZtBgwZ1qoYPPvhA1dXVrcZaOhwOLV26VGvWrNHnn3/e4WMROgAAMJvR2Sde7h8fH6/4+HiDJ/XstttuU1paWqt1N9xwg2677TbNnTvXq2MROgAAMJtJ3Stdoby8XDU1NSovL5fD4VBxcbEkafDgwerRo4ckKSUlRdnZ2ZoxY4Z69eqlXr16tTrGeeedp4SEBA0ZMsSrcxM6AAA4hyxbtkwvvfSS+/Vll10mSXr//fc1ceJESVJpaanq6upMPzehAwAAs/lxS0dubu5Zb0txtpml3ozj+D5CBwAAZmuW5DSwv5F9/Ri3QQcAAD5BSwcAAGZzyFj3SpC2dBA6AAAwW7OM9SUQOgAAQIcQOjxiTAcAAPAJWjoAADBbk2jp8IDQAQCA2ZwyNpA0SB/ATvcKAADwCVo6AAAwW7MkI0+BD9KWDkIHAABmI3R4RPcKAADwCVo6AAAwW5No6fCA0AEAgNkcInR4QPcKAADwCVo6AADoCkHaWmEELR0AAMAnCB0AAMAnCB0AAMAnCB0AAMAnGEgKAIDpmr5bjOwffGjpAAAAPkFLBwAApmv+bjGyf/AhdAAAYDq6VzyhewUAAPgELR0AAJiO7hVPCB0AAJiuWca6SIIzdNC9AgAAfIKWDgAATMdAUk8IHQAAmI4xHZ4QOgAAMB1jOjxhTAcAAPAJWjoAADAd3SueEDoAADAdA0k9oXsFAAD4BC0dAACYju4VTwgdAACYjtkrntC9AgAAfILQAQCA6ZpNWLrGypUrNX78eHXv3l2xsbEd3q+kpEQ//elPFRMTo+joaF1++eUqLy/36tyEDgAATNdkwtI1Tp06pZkzZ2rhwoUd3ufw4cOaMGGCUlJStG3bNv3jH//Qgw8+qKioKK/OzZgOAADOIStWrJAk5ebmdnif3/zmN0pPT9fjjz/uXnfBBRd4fW5aOgAAMJ053St2u73V0tjY6OPPITmdTr399tu66KKLdMMNN6hPnz5KTU3Vxo0bvT4WoQMAANO1zF7p7HI6dCQlJSkmJsa9ZGdn+/hzSNXV1aqvr9djjz2myZMn67//+781Y8YM/eu//qu2b9/u1bHoXgEAwHTm3KejoqJCNpvNvTYyMtLj1llZWVq1atUZj1hSUqKUlBSvK3E6nZKkadOm6Z577pEkjRw5Uh999JHWrVunq6++usPHInQAAOCnbDZbq9DRnqVLl2rOnDln3GbQoEGdqqF3794KDw/XsGHDWq0fOnSoPvzwQ6+ORegAAMB0vn32Snx8vOLj4w2cr30RERG6/PLLVVpa2mr9Z599pgEDBnh1LEIHAACm898HvpWXl6umpkbl5eVyOBwqLi6WJA0ePFg9evSQJKWkpCg7O1szZsyQJN1777265ZZbdNVVV+maa67R5s2b9de//lXbtm3z6tyEDgAAziHLli3TSy+95H592WWXSZLef/99TZw4UZJUWlqquro69zYzZszQunXrlJ2drbvuuktDhgzRm2++qQkTJnh1bkIHAACm898HvuXm5p71Hh0ul6vNup///Of6+c9/bujchA4AAEzHA9884T4dAADAJ2jpAADAdP7bvWIlQgcAAKZrkrGf2K6bvWIlulcAAIBP0NIBAIDp6F7xhNABAIDpmL3iiV90r6xdu1bJycmKiopSamqqCgoKrC4JAAADzHm0fbCxPHS89tprWrJkiR566CHt27dPI0aM0A033KDq6mqrSwMAACayPHT8/ve/1/z58zV37lwNGzZM69atU/fu3fXHP/7R6tIAAOikJhOW4GNp6Dh16pSKioqUlpbmXhcaGqq0tDTt2rXLwsoAADCC7hVPLB1I+uWXX8rhcKhv376t1vft21cHDhxos31jY6MaGxvdr1seRlPetWUGraPf+3vKykIC0PHv/nLtvMe1M6bl+n0hqfFMG6KNlt8KT88VMZ/R/3WC83/dgJq9kp2drRUrVrRZ/5oFtQSTv1pdQADj2nUe186Yl60uIIB99dVXiomJ6ZJjR0REKCEhQVVVqw0fKyEhQRERESZU5T8sDR29e/dWWFiYjh071mr9sWPHlJCQ0Gb7Bx54QEuWLHG/rq2t1YABA1ReXt5lX6BgZrfblZSUpIqKCtlsNqvLCShcu87j2hnD9eu8uro69e/fX3FxcV12jqioKJWVlenUKePteBEREYqKijKhKv9haeiIiIjQ6NGjtWXLFk2fPl2S5HQ6tWXLFmVmZrbZPjIyUpGRkW3Wx8TE8C+fATabjevXSVy7zuPaGcP167zQ0K4dzhgVFRV0YcEslnevLFmyRLNnz9aYMWM0duxYrVmzRg0NDZo7d67VpQEAABNZHjpuueUWHT9+XMuWLVNVVZVGjhypzZs3txlcCgAAApvloUOSMjMzPXannE1kZKQeeughj10uODuuX+dx7TqPa2cM16/zuHbWC3H5Zu4QAAA4x1l+R1IAAHBuIHQAAACfIHQAAACfIHQAAACfCOjQsXbtWiUnJysqKkqpqakqKCiwuqSAsGPHDk2dOlWJiYkKCQnRxo0brS4pIGRnZ+vyyy9Xz5491adPH02fPl2lpaVWlxUwcnJyNHz4cPdNrcaNG6e//e1vVpcVkB577DGFhIRo8eLFVpcSEJYvX66QkJBWS0pKitVlnZMCNnS89tprWrJkiR566CHt27dPI0aM0A033KDq6mqrS/N7DQ0NGjFihNauXWt1KQFl+/btWrRokXbv3q38/Hw1NTXp+uuvV0NDg9WlBYTzzz9fjz32mIqKirR3715de+21mjZtmj799FOrSwsohYWFeu655zR8+HCrSwkoF198sY4ePepePvzwQ6tLOicF7JTZ1NRUXX755Xr22Wclnb59elJSku68805lZWVZXF3gCAkJ0YYNG9y3oUfHHT9+XH369NH27dt11VVXWV1OQIqLi9N//Md/aN68eVaXEhDq6+s1atQo/eEPf9Cjjz6qkSNHas2aNVaX5feWL1+ujRs3qri42OpSznkB2dJx6tQpFRUVKS0tzb0uNDRUaWlp2rVrl4WV4VxSV1cnSV368Khg5XA49Oqrr6qhoUHjxo2zupyAsWjRIt14442t/r8PHXPw4EElJiZq0KBBysjIUHl5+dl3gun84o6k3vryyy/lcDja3Cq9b9++OnDggEVV4VzidDq1ePFiXXHFFbrkkkusLidg7N+/X+PGjdPJkyfVo0cPbdiwQcOGDbO6rIDw6quvat++fSosLLS6lICTmpqq3NxcDRkyREePHtWKFSt05ZVX6pNPPlHPnj2tLu+cEpChA7DaokWL9Mknn9Av7KUhQ4aouLhYdXV1Wr9+vWbPnq3t27cTPM6ioqJCd999t/Lz83l6aSdMmTLF/c/Dhw9XamqqBgwYoNdff52uPR8LyNDRu3dvhYWF6dixY63WHzt2TAkJCRZVhXNFZmamNm3apB07duj888+3upyAEhERocGDB0uSRo8ercLCQj311FN67rnnLK7MvxUVFam6ulqjRo1yr3M4HNqxY4eeffZZNTY2KiwszMIKA0tsbKwuuugiHTp0yOpSzjkBOaYjIiJCo0eP1pYtW9zrnE6ntmzZQv8wuozL5VJmZqY2bNigrVu3auDAgVaXFPCcTqcaGxutLsPvTZo0Sfv371dxcbF7GTNmjDIyMlRcXEzg8FJ9fb0OHz6sfv36WV3KOScgWzokacmSJZo9e7bGjBmjsWPHas2aNWpoaNDcuXOtLs3v1dfXt0r4ZWVlKi4uVlxcnPr3729hZf5t0aJFysvL01/+8hf17NlTVVVVkqSYmBh169bN4ur83wMPPKApU6aof//++uabb5SXl6dt27bp3Xfftbo0v9ezZ882Y4eio6PVq1cvxhR1wK9//WtNnTpVAwYM0JEjR/TQQw8pLCxMs2bNsrq0c07Aho5bbrlFx48f17Jly1RVVaWRI0dq8+bNbQaXoq29e/fqmmuucb9esmSJJGn27NnKzc21qCr/l5OTI0maOHFiq/Uvvvii5syZ4/uCAkx1dbVuv/12HT16VDExMRo+fLjeffddXXfddVaXhiD3z3/+U7NmzdJXX32l+Ph4TZgwQbt371Z8fLzVpZ1zAvY+HQAAILAE5JgOAAAQeAgdAADAJwgdAADAJwgdAADAJwgdAADAJwgdAADAJwgdAADAJwgdAADAJwgdAADAJwgdAADAJwgdQJA4fvy4EhIS9Lvf/c697qOPPlJERESrJzIDgFV49goQRN555x1Nnz5dH330kYYMGaKRI0dq2rRp+v3vf291aQBA6ACCzaJFi/Tee+9pzJgx2r9/vwoLCxUZGWl1WQBA6ACCzbfffqtLLrlEFRUVKioq0qWXXmp1SQAgiTEdQNA5fPiwjhw5IqfTqc8//9zqcgDAjZYOIIicOnVKY8eO1ciRIzVkyBCtWbNG+/fvV58+fawuDQAIHUAwuffee7V+/Xp9/PHH6tGjh66++mrFxMRo06ZNVpcGAHSvAMFi27ZtWrNmjV5++WXZbDaFhobq5Zdf1gcffKCcnByrywMAWjoAAIBv0NIBAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB8gtABAAB84v8D2sbdDyuMeekAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pp.plot_grid(g, cell_value=u[1::2], plot_2d=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand the inner workings of the discretization, and to recover the traction on the grid faces, some more detail is needed. The mpsa discretization creates two sparse matrices \"stress\" and \"bound_stress\". They define the discretization of the cell-face traction:\n",
    "\\begin{equation}\n",
    "T = \\text{stress} \\cdot u + \\text{bound\\_stress} \\cdot u_b\n",
    "\\end{equation}\n",
    "Here, $u$ is a vector of cell center displacement and has length g.dim $*$ g.num_cells. The vector $u_b$ is the boundary condition values. It is the displacement for Dirichlet boundaries and traction for Neumann boundaries and has length g.dim $*$ g.num_faces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The global linear system is now formed by momentuum balance on all cells. A row in the discretized system reads\n",
    "\\begin{equation}\n",
    "-\\int_{\\Omega_k} f dv = \\int_{\\partial\\Omega_k} T(n)dA = [div \\cdot \\text{stress} \\cdot u + div\\cdot\\text{bound\\_stress}\\cdot u_b]_k,\n",
    "\\end{equation}\n",
    "\n",
    "The call to mpsa_class.assemble_matrix_rhs(), creates the left-hand side matrix $ \\text{div} \\cdot \\text{stress} $, the right-hand side vector, which consists of $\\text{f}$ and $-\\text{div} \\cdot \\text{bound\\_stress}$ (note sign change).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also retrieve the traction on the faces, by first accessing the discretization matrices `stress` and `bound_stress`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGzCAYAAAAIWpzfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABrAUlEQVR4nO3deVxU5f4H8M+wiwIuyKaAuYA7mgvilnoxUn+Wtmhq7lEpVIrmmltlWJq3zbQsl256NW9pXTW9aqGpkCsqbrghKouCLMOwM+f3h6+ZHEWZgTNz5oHP+/Wal8w485wPM2eY7zzPc56jkiRJAhEREZEAbJQOQERERGQsFi5EREQkDBYuREREJAwWLkRERCQMFi5EREQkDBYuREREJAwWLkRERCQMFi5EREQkDBYuREREJAwWLkREAJKSkqBSqbBu3bpKPX7dunVQqVRISkqSNRcRGWLhQiSoK1eu4PXXX0fTpk3h5OQEV1dX9OjRA5999hkKCgqUjqeXkpKChQsXIj4+3qTHXbt2DZGRkQgICICzszOcnZ3RunVrRERE4PTp0+YJS0RWz07pAERkuh07duCll16Co6MjxowZg7Zt26K4uBgHDx7EO++8g7Nnz+Kbb75ROiaAe4XLokWL0KRJE3To0MGox2zfvh3Dhw+HnZ0dRo0ahaCgINjY2ODChQv4+eefsXLlSly7dg3+/v6y5fT390dBQQHs7e1la5OI5MfChUgw165dw8svvwx/f3/8/vvv8Pb21v9fREQELl++jB07dlg0k0ajQe3atWVp68qVK/rfb9++fQa/HwB89NFH+Oqrr2BjY3qH8eNyqlQqODk5VSozEVkOh4qIBPPxxx8jLy8P33333UMf6gDQvHlzvP32249t488//8RLL70EPz8/ODo6wtfXF1OnTjVqiEk3l2P//v2YPHkyPDw80Lhx43LvGxMTgy5dugAAxo8fD5VKVeE8ko8//hgajQZr164t9/ezs7PDW2+9BV9fX9lyAqbNcTl79iz69euHWrVqoXHjxvjggw+g1WorfBwRVR17XIgE89///hdNmzZF9+7dK93Gli1bkJ+fj0mTJqFBgwY4cuQIvvjiC9y8eRNbtmwxqo3JkyejYcOGmD9/PjQaTbn3adWqFd577z3Mnz8fr732Gnr16gUAj82+fft2NG/eHMHBwab/YpXMaYq0tDT07dsXpaWlmDVrFmrXro1vvvkGtWrVkiEtEVWEhQuRQHJzc3Hr1i0899xzVWrno48+Mvigfe2119C8eXPMmTMHycnJ8PPzq7CN+vXrY9++fbC1tX3kfTw9PTFgwADMnz8fISEheOWVVx7bZm5uLlJSUjBkyJCH/i87OxulpaX667Vr1zaqWDAmpyk++ugj3LlzB3/99Re6du0KABg7dixatGghS/tE9HgcKiISSG5uLgDAxcWlSu3c/4Gv0WiQkZGB7t27Q5IknDx50qg2wsPDZSsGdHS/X506dR76vz59+qBhw4b6y4oVKxTJuXPnTnTr1k1ftABAw4YNMWrUKNm2QUSPxh4XIoG4uroCANRqdZXaSU5Oxvz58/Hrr78iKyvL4P9ycnKMauOJJ56oUoby6AqyvLy8h/7v66+/hlqtRnp6eoU9N/eTO+f169fLHcYKDAyUdTtEVD4WLkQCcXV1hY+PDxISEirdRllZGfr374+7d+9i5syZaNmyJWrXro1bt25h3LhxRk8yNcecDjc3N3h7e5f7++mKBVMXeOPcE6LqhUNFRIL5v//7P1y5cgWxsbGVevyZM2eQmJiITz75BDNnzsRzzz2H0NBQ+Pj4yJz0HpVKZdL9Bw0ahMuXL+PIkSNmyVNV/v7+uHTp0kO3X7x4UYE0RDUPCxciwcyYMQO1a9fGq6++ivT09If+/8qVK/jss88e+XjdfA9JkvS3SZL02MdUhW7dlOzsbKPuP2PGDDg7O2PChAnl/n7351bCwIEDERcXZ1BY3blzBxs2bFAwFVHNwaEiIsE0a9YMGzduxPDhw9GqVSuDlXMPHz6MLVu2YNy4cY98fMuWLdGsWTNMnz4dt27dgqurK3766aeH5rrImbdu3bpYtWoVXFxcULt2bQQHBz9y7kmLFi2wceNGjBgxAoGBgfqVcyVJwrVr17Bx40bY2Ng8dk0Wc5oxYwb+9a9/4ZlnnsHbb7+tPxza39+fpyIgsgSJiISUmJgohYeHS02aNJEcHBwkFxcXqUePHtIXX3whFRYWPvax586dk0JDQ6U6depI7u7uUnh4uHTq1CkJgLR27drHPnbt2rUSAOno0aNGZ/3ll1+k1q1bS3Z2dkZtQ5Ik6fLly9KkSZOk5s2bS05OTlKtWrWkli1bSm+88YYUHx9f4eNNzXnt2jWjs50+fVp66qmnJCcnJ6lRo0bS+++/L3333XcSAOnatWtGbY+IKkclSQr3uxIREREZiXNciIiISBgsXIiIiEgYLFyIiIhIGCYXLgcOHMDgwYPh4+MDlUqFbdu2VfiYmJgYPPnkk3B0dETz5s2NOvsqERER0YNMLlw0Gg2CgoKMPk/ItWvXMGjQIPTt2xfx8fGYMmUKXn31VezevdvksERERFSzVemoIpVKha1bt5Z7JledmTNnYseOHQZLeL/88svIzs7Grl27KrtpIiIiqoHMvgBdbGwsQkNDDW4LCwvDlClTHvmYoqIiFBUV6a9rtVrcvXsXDRo0MHn5cCIiIlKGJElQq9Xw8fGBjY0802rNXrikpaXB09PT4DZPT0/k5uaioKCg3BOgRUdHY9GiReaOJgwutUNEJLaa/qX7xo0bsq12bZVL/s+ePRtRUVH66zk5OfDz88ONGzfg6uqqYDIiIiLT5eTkKB1BEbm5ufD19YWLi4tsbZq9cPHy8nroRGnp6elwdXV95OnmHR0d4ejo+NDtrq6uLFyIiEg4Nf2zS84eJ7Ov4xISEoJ9+/YZ3LZnzx6EhISYe9NERERUzZhcuOTl5SE+Ph7x8fEA7h3uHB8fj+TkZAD3hnnGjBmjv/8bb7yBq1evYsaMGbhw4QK++uor/Pjjj5g6dao8vwERERHVGCYXLseOHUPHjh3RsWNHAEBUVBQ6duyI+fPnAwBSU1P1RQwAPPHEE9ixYwf27NmDoKAgfPLJJ/j2228RFhYm069ARERENYUQZ4fOzc2Fm5sbcnJyavw4IRERkSjM8fnNcxURERGRMFi4EBERkTBYuBAREZEwWLgQERGRMFi4EBERkTBYuBAREZEwWLgQERGRMFi4EAAgKSkJKpXqoUufPn2UjkZERKRnlWeHJsvz9fVFamqq/npaWhpCQ0PRu3dvBVMREREZ4sq59JDCwkL06dMHDRs2xC+//AIbG3bMERGR6czx+c0eF3rIhAkToFarsWfPHhYtRERkVVi4kIEPPvgAu3fvxpEjR+Di4qJ0HCIiIgMsXEjvp59+wnvvvYfffvsNzZo1UzoOERHRQ1i4EAAgISEBY8aMwcyZM9GmTRukpaUBABwcHFC/fn2F0xEREd3DCQwEADh27Bjy8/PxwQcfwNvbW395/vnnlY5GRESkx6OKiIiIyCzM8fnNHhciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBgsXIiIiEgYLFyIiIhIGCxcy2ooVK9CkSRM4OTkhODgYR44ceez9t2zZgpYtW8LJyQnt2rXDzp07H3nfN954AyqVCp9++qnMqYmIqDph4UJG2bx5M6KiorBgwQKcOHECQUFBCAsLw+3bt8u9/+HDhzFixAhMnDgRJ0+exJAhQzBkyBAkJCQ8dN+tW7ciLi4OPj4+5v41iIhIcCxcyCjLly9HeHg4xo8fjz/++AO7d+9GdnY2OnXqVG7Py2effYZnnnkG77zzDhISErBlyxaUlpaid+/eBj0vSUlJGDNmDDQaDVJTU7FgwQKMGTMGKSkplvz1iIhIECxcqELFxcU4fvw4QkND9T0vCxcuxLPPPgtbW9tye15iY2MRGhpq0PPyxhtvwN7eXt/zotVqMWbMGPj4+GDx4sXw9vbG+PHjcfHiRTz77LMK/bZERGTNWLhQhTIyMlBWVgZPT0+DnpfAwEAAgFqtho+Pj8G8l7S0NHh6ehr0vLRt2xZFRUWwtbVFUFAQfHx8kJeXhwsXLmDYsGGwt7eHv78/2rVrh+PHj8PJyQmhoaG4dOmSkr8+ERFZERYuZLSSkhJ9zwsAXLhwAcnJyQgODkbv3r3Lnfei63kBgKtXryInJwc9e/aEr68v8vLycObMGZw9e1Z//99//x2bNm2CSqXC77//jtq1ayMsLAyFhYWW/WWJiMgqsXChCrm7u8PW1haXL1/W97wA9ybg+vv7o0ePHtBoNFi1ahWcnZ2xZs0aeHl5IT09Xd/zAgC//fYbXFxcMHToUGRmZiI/Px+lpaVo37497OzscP36dWzfvh2FhYUYMWIEunfvju+//x4pKSnYtm2bgs8AERFZCxYuVCEHBwd06tQJBw8e1N9WWFiIzMxM9OzZU3+bjY0NQkNDERsbi06dOmHp0qUoKSnBxIkTMXHiRFy+fBlt27YFADg5OeH06dN444030KJFC8THx6Nhw4YAAK1Wi23btuGFF15AYWEhgoODERsbi1OnTmHEiBHw9fVFrVq10KpVK3z22WeWfTKIiEhRLFzIKFFRUdi4cSNUKhWOHTuGV199FQDw8ssvIz09HSkpKZg9ezY8PT2RlpaG1NRUpKSkoE6dOnj++eexdetWFBcXY8iQIUhPT4ePjw/atm2Ldu3aITs7G4GBgbh79y4A4Ntvv8WBAweQkpKC559/Xt/m8ePH4eHhgR9++AFnz57F3LlzMXv2bHz55ZdKPjVERGRBdkoHIDEMHz4cd+7cwbRp0zBlyhQEBQUBAOrWrYt9+/bBwcEBqamp8PDwQEFBAc6cOYPo6GgsXrwYGzZsgEqlAgBs27YNZWVlCAkJ0bctSRKGDh2KsrIyAEDfvn2xbNkynD9/Hjk5OWjYsCFCQkIwYcIEg0xNmzbFH3/8gWnTpuHNN99EVlYW6tata5knhIiIFMEeFzJaZGQkvv/+e9jb2yMiIgI2NjZ47733oNFoEBsbi3Xr1uGnn35Ceno66tati1mzZqFNmzYAgMmTJwMAjh8/jmPHjiEyMhIAkJKSgqKiIoO1YN566y1s27YNmzdvhpeXF3Jzc3H06NFyM/3vf/9DvXr1zPybExGRtWDhQiYZPnw4li1bhvfeew+SJOHUqVPYtWsXPD09odVqcfPmTdjZ2cHDwwPnz5/HX3/9hejoaOzZswfAvd4VrVaL+vXrAwB27tyJ3Nxc3LlzR7+N7du3IyUlBY6OjvDx8UFpaSlSU1MRFxdnkGX69Om4desWoqKiLPcEEBGRoli4kNFmzZoFlUqFN998E8nJyZAkCenp6ejWrRtUKhVsbW3h5OSEsWPHAgAmTpwIR0dHzJo1CxcvXkTdunVRVlYGlUqFn376CQsXLkRCQgLOnDmDDRs2wMHBQT93ZsOGDWjQoAGuXr2KOnXqwNfXF7Gxsfosv/zyC5YvX46oqCh07dpVkeeDiIgsj4ULGW3atGk4f/68weXdd9+Ft7c37O3t0b59e+zYsQPNmzfH7du3kZaWBkdHRwBAaWkp1Go1oqKiYGNjg6ioKPznP//Btm3b0LZtW3h5eaG4uBhdunTR/3+XLl1QUFCASZMmwcvLC2lpaQCA+Ph4vPjiixg8eDCWLVum5FNCREQWxsKFjNawYUO0bNkS69atQ6tWrdCqVSt88MEHSE1NRUlJCU6fPo1evXph4sSJyM7OxoABA+Dl5QXg3sJyWq0WU6dORf369fH5558jISEBAwcOBAB06tQJ9vb2OHfuHOzs7JCWloZTp06hqKgIzz33nD7D2bNn0aNHDzRv3hy//PKLIs8DEREph0cVkcmmTZuGcePGPfY+b7/9NrZv3467d+/i0KFDiIyMxMsvvwwPDw/cvXsX9vb2aNmyJb7//nt07doVbm5umDhxIv7973+juLgYMTExmDNnDkJCQtCtWzekp6dDkiT07dsXjo6OuHjxImxtbQHcmzcD3Fsob+7cuVi0aJG5nwIiIlIICxcyWcOGDfWLxT3Kv//9b4wZMwY7duxAWFgYhg0bhs8//1zf89KuXTtcvHgR+fn5+sf885//RGlpKb799ls888wzGDhwIL766ivs3r0bycnJuH37tsEkXl3BUr9+fdy9exd79+7FE088gYKCAvP84jIrKirSD6WJRpIkSJIEGxsxO221Wq2w2cvKyiBJEuzsxPzzXVBQAEdHR4PnX5Ik2Nvbw97eXsFkJAox93yyevXr18f27dsxYMAApKen44033sCpU6f0PS/BwcGQJAm3bt0y6HlZvXo17OzssHPnTkRGRuLixYt45513YGdnh9q1a+PKlSuQJAllZWW4dOkSfHx8kJSUhOeffx6NGjVCrVq1oFarlf71K5STk4NLly6hY8eO+p4jkVy+fBlOTk5o3Lix0lFMVlpaiqNHj6Jz585CflAmJycjPz8fLVu2VDqKyQoLC5GQkAB/f3/9lx9JkqBSqVCrVi0hXw+yPJWk+9pqxXJzc+Hm5oacnBy4uroqHYdMcPfuXURGRuK///0vbGxs8MILL+Dzzz9HnTp1AABJSUl44okn8Mcff6BPnz4A7v1xmzZtGjZu3Ijs7OxHtr1x40Y4OzvjzJkzmDdvHn744Qd9u0QkDhsbG/Tr1w+1a9dWOgrJzByf3yxcyGpJkgS1Wq3/9/Lly3jppZcwdOhQREZGonbt2tBoNLh27RqaNm2qnwgsiry8PJw+fRpPPvkknJyclI5jkqysLFy5cgWdO3dWOorJUlJSkJ2djdatWysdxWQXLlyAra0tWrRooXQUkxQVFSEhIQF169ZF06ZNoVKpUFpaCkmSYGtri7i4OISGhrJwqYZYuLBwqZF0r/+j6HpeiEgsNjY26NWrF/bv38/CpZoyx+c357iQ1XNxcUFOTo7+emJiIgYNGoSXXnoJ06dPh6urK7Kzs3H27FkEBgZWOHHYmhQVFeHo0aNo27atcOdZiouLQ8uWLYXLffz4cfj7+8Pd3V3pKEaTJAknTpyAh4cHfH19lY5jNN0yCbqzuet6WnQLUcbFxenPY6bVahVOS6IQc1o91SgqlQqurq5wcXEBAAQEBGDHjh3YsmULPvnkE+Tm5sLW1hYtW7bEhQsXkJmZqXBi4zk6OsLf318/6VgkderUQV5entIxTKLVaqHRaIT7Zp+amorS0lI0atRI6ShGKy0tRUJCAhwdHdGyZUuDosXBwUHfS6r7UpKRkaFkXBIIh4pIGBwyIqoebGxs8NRTT+GPP/6ASqWCJEkcKqqmOFRENZpuyEg3WRe4N2w0bNgwDB06VH+6AN2E3TZt2ghT6GZmZuLy5ct48sknhTkkNDMzE8nJyejYsaPSUYyWlpaGjIwMtG3bVukoRrt27Rry8/PRunVr/bCKNSsrK8OFCxcgSRJatWoFW1tbg54We3t7HDhwQN8z2qJFCyQmJiqcmkTCHhcSDnteiMRUXk9Lnz59EBMTwx6Xaoo9LkR4fM/LkCFD4O7ujvz8fOTn5yM9PR3Nmze3+sONs7KyUKtWLSQmJiIwMFCIFXUlSUJsbCw6dOggTKF4+fJl1K9fHyqVCvXq1VM6ToWuXbsGe3t7NG7cGFlZWVadWavVIjk5GcXFxWjatCns7OxQVlYGZ2dnBAYG6ntadCdLbd26Nc6ePatwahIRCxcSjm6ybm5u7kNHWHz33Xf47rvvFEpGROVZtWoVRowYAeDeWjQAhDr6j6wLCxcS1oOHSQPApk2b8PrrryM6OhpdunRBSUkJzp07hw4dOigT0gh79+5FdHQ0Zs+ejaCgIKjVajRt2lTpWBXS5Y6MjMTQoUOVjmOU+Ph4pKWlYenSpZg9ezZCQ0OVjvRIarUaN27cQOvWrQ32EWvNfOvWLQAwOPLp4sWLmDx5MjQaDW7cuAEAaN++PU6dOqVIRqoeWLiQsHQ9L/fTnb+lY8eO6NatGwoKClBUVITg4OAqbSs3Nxf/+te/MHnyZNknSCYlJQEAmjdvjg4dOiAjI0O2Quuvv/5CXl4e/vGPf8jS3v10uf38/Kr8/D5o/fr16N+/P3x8fGRtNy8vTz+s1bx5c1lz5+Tk4IcffpBtH8nMzISzszOCg4MN9hG5n+u//voLGo0G/fr1q1I7ugm2AQEBAKBfFVfn6tWrAGDVw10kBq7jQtWSSqXSXx68XpnLn3/+iZkzZyIzM7PKbT14qSh/VS7ffPMNPvnkE9kz35/bHG1PmTIFO3fuNPtzLWfbun0kKytL1rwP5pb7OVm1ahWWLVtmlv3hftbc60liYY8LCe3+CbrAvZ4RALhx4wYSExPh4OAA4N4CXlVx9+5dAEB6ejpKS0ur1NaDdCeSzM7ORk5ODgoLC6ucV0fX4yRXe/fT5c7PzzdL+zk5ObK3e//+kp2dLWv7WVlZAO7tI8XFxVVuLy8vD6WlpUhNTTXYR+R+TgoLC1FcXFzldjUaDYC/32uSJCE5OVn//zwilORSqcJlxYoVWLp0KdLS0hAUFIQvvvgCXbt2feT9P/30U6xcuRLJyclwd3fHiy++iOjoaKs/0oOsn1qtLvfQ6FdffRUA8N///hcAcPPmzSptR7fmRGpqKgoKCqrU1oN0RdHdu3eRnZ2NkpKSKufVyc/PR1FRkWzt3U+XOz8/X/b2JUlCdna2WXLrioC7d+/K2r5u5deUlBT9h3hVlJaWorS0FDdv3jTYR+R+TuTaR3Tvi/vbKSoqqlKbROUxuXDZvHkzoqKisGrVKgQHB+PTTz9FWFgYLl68CA8Pj4fuv3HjRsyaNQtr1qxB9+7dkZiYiHHjxkGlUmH58uWy/BJUcz04QffQoUMYOHAgtm7dim7dusHOzg6HDh1Cly5dqrSd9PR0APe6u+U+GuLKlSsAgKZNm8Lf3x8ZGRmyLerWoEEDFBcXV/n3L48ut7u7u+ztq1Qq+Pv7y97u3r179UeiNW3aVNb2dYf5duzYEQ0aNKhye5mZmbhw4QK6dOlisI/I/ZzUr18fJSUlVW734sWLAIDAwEAAEO4UFiQOk+e4LF++HOHh4Rg/fjxat26NVatWwdnZGWvWrCn3/ocPH0aPHj0wcuRINGnSBE8//TRGjBiBI0eOVDk8kUp1b4Ku7qJbwMrZ2Rmurq7lzmsgIiJxmVS4FBcX4/jx4waH49nY2CA0NBSxsbHlPqZ79+44fvy4vlC5evUqdu7ciYEDBz5yO0VFRcjNzTW4UPW0cOHChyb56Y4MIiIiepBJQ0UZGRkoKyuDp6enwe2enp76RYUeNHLkSGRkZKBnz56QJAmlpaV44403MGfOnEduJzo6GosWLTIlGgmsTZs22Lt3r/66nR3njBMRUfnMfjh0TEwMPvzwQ3z11Vc4ceIEfv75Z+zYsQPvv//+Ix8ze/Zs5OTk6C+6hYuoerKzs4OXl5f+4u7urnQkIiKyUiZ9tXV3d4etra1+oqJOeno6vLy8yn3MvHnzMHr0aP1RHu3atYNGo8Frr72GuXPnwsbm4drJ0dFRiHO1kDwuXboEHx8fODk5ISQkBNHR0fDz81M6FhERWSGTelwcHBzQqVMn7Nu3T3+bVqvFvn37EBISUu5j8vPzHypOdKspctY5BQcHY926ddi1axdWrlyJa9euoVevXgZrsxAREemYPJkgKioKY8eORefOndG1a1d8+umn0Gg0GD9+PABgzJgxaNSoEaKjowEAgwcPxvLly9GxY0cEBwfj8uXLmDdvHgYPHmywHDTVTAMGDND/3L59ewQHB8Pf3x8//vgjJk6cqGAyIiKyRiYXLsOHD8edO3cwf/58pKWloUOHDti1a5d+wm5ycrJBD8u7774LlUqFd999F7du3ULDhg0xePBgLF68WL7fgqqNunXrIiAgAJcvX1Y6ChERWaFKHb4RGRmJyMjIcv8vJibGcAN2dliwYAEWLFhQmU1RDZOXl4crV65g9OjRSkchIiIrxJMskqKmT5+O/fv3IykpCYcPH8bQoUNha2uLESNGKB2NiIisEBfMIEXdvHkTI0aMQGZmJho2bIiePXsiLi5O9mX1iYioemDhQoratGmT0hGIrJZGo8HVq1eVjkFkVThURERkhbKzs/Hmm2+iZ8+eKCkpUToOkdVg4UJE1dqhQ4fw448/CrVulEajwSuvvIKUlBTk5OQgPj5e6UhEVoOFCxEZJTMzE3fu3FE6hskWLFiA8ePHY/To0cjPz1c6jlH++usvnD59Wn/9t99+UzANkXXhHBcihd26dUvpCEYZMWIEDh8+jA4dOqBr165ISUlROlKFJEnCmTNnAABbt27Fb7/9hh49eiicqmL9+vXD9u3b8eOPP+Lw4cNc14joPixciBSg+0DdsGEDVq1apXQco3h7ewMATp48iZMnT8LBwUHhRBVLT09HXl4egHvnQHNxcRFmsmv79u3h4OCAlStX4rXXXsOGDRuUjkRkFThURGRhCQkJaNasGUJCQvDll18KUQAAQNu2bSFJEmxsbPDee+/hyy+/VDpShTw8PLB06VJs3boVKSkp+PbbbzFjxgylY5msX79+SkcgshosXIgszMXFBfn5+VCpVFCpVBg0aJDSkYzyzDPPoEuXLti9ezemTZsmxLnGbGxsMHnyZDz99NNwcnJSOg4RyYCFC5GF+fv749dff4WDgwNUKhW6deumdCSjtGvXDjExMejevbvSUYioBuMcFyIFdO3aFTt27EBiYiIcHR2VjkNEJAz2uBApJCQkBGPHjlU6BhGRUFi4EBERkTBYuBAREZEwWLgQERGRMFi4EBERkTBYuBAREZEwWLgQERGRMFi4EBERkTBYuBAREZEwWLgQERGRMFi4EFVg7dq1iIqKAgCEhYXhr7/+UjiRcV5++WXs3LkTR48exYABA1BcXKx0pApdvHgR3bt3R0lJCaKjo/HRRx8pHcko3333nX4f6d+/P44cOaJwIuMMHz4cu3bt0u8jJSUlSkciqhDPVURUgaysLNy6dQvAvQ/WsrIyhRMZ5/r161Cr1QCAxMRE2NhY//cUOzs7nD59GgCQnp6OjIwMhRMZJzs7GykpKQDu7SNarVbhRMZJSkrS7yOXLl2CSqVSOBFRxaz/LxmRwsLDw+Hs7AwACA4OFubsyO+++67+5zlz5sDOzvq/pzRr1gwvvfQSAMDW1hZTp05VOJFxXnvtNdSqVQvAvXNQiXLG73nz5ul/FmUfIWLhQlQBFxcXvPDCCwCAGTNmKJzGeAMHDkS9evXg4OCA0aNHKx3HaLNnzwYA9OzZEz4+PgqnMY6o+8igQYNQt25dODo64pVXXlE6DpFRWF4TGeHzzz9H79698cwzzygdxWgqlQq7du2CWq2Gg4OD0nGMFhAQgH/961/o27ev0lFM8sUXX6BPnz54+umnlY5iNJVKhd27dyMvL0+xfUQ3rFZUVITatWsrkoHEwsKFyAgODg4YOXKk0jFM1rZtW6UjVMrzzz+vdASTOTg4YMSIEUrHMJlS+4gkSQCgn8gswuRxsg4cKiIiIos7deoUAMDX1xfAveE2ImOwx4Wqpd9//x03b95EcXExEhMTkZmZqXSkR4qLi9P/e/fuXajValy8eFHhVBXT5T5x4gQ2bdqkcBrjnD17FllZWQD+zm+t8vLykJKSghs3bhjsI9YqNTUVAHDy5En9bVevXn3ofomJiQCAOnXqIDs7G40aNcKlS5csE5KqBZWk66+zYrm5uXBzc0NOTg5cXV2VjkNW7JdffsGQIUOUjkFE91m4cCGef/55XL16FfXq1UNWVhb69OmDmJgY/b+hoaGc41INmePzmz0uVK3Uq1cPADBz5kwEBAToe1ysea5HXFwcVq9ejfDwcAQEBECtVsPf31/pWBXS5R4xYgRCQ0OVjmMUXY/L2rVrER4ebtWHLet6XAICAgz2EWvNrOtx8fb21t929epVLF68GC4uLvr/DwoKQkxMjBIRqZpg4UJCkyRJv4AWcG8hMADw8vJCkyZN4OjoCG9vb3h5eVVpO8nJyVi0aBG++uorODo6VqmtB924cQMA0KhRI3Tq1AlFRUX6AqyqtmzZgrS0NLz55puytHc/Xe4nnngCbdq0kbXtKVOmYOzYsejYsaOs7TZo0EA/lNGoUSNZc1++fBkfffQRVq5cKct6KEVFRWjatCkaNmxosI/I/Vxv2LABGo0Gr732WpXa0c1Vuf9bdWlpqf7nLl264MCBA1XaBhHAwoUEp1ar4ebm9tDtuoXL/vzzTwD3xtOr4vbt2zh79iy0Wm2V23qQk5OT/l9HR0eUlpbKto0rV64gKSlJ9szA37nt7e1lb//YsWMIDQ1Fr169ZG0XgP6wXycnJ1lzp6am4syZM7C3t5dlyMPGxgZ5eXmoU6eOwT4i93N98eJFqNXqKrdbUFAAwPC9VrduXf3PIqzcTGJg4UJCc3FxQU5Ojv76oUOHMHDgQGzduhXdunWDnZ0dDh06hCeeeKJK29FNlvX19TXoCpdDw4YN9f82aNAAkiRVOa+Oi4sLHBwcZGvvfrrcderUkb19lUoFd3d32du9cuUK6tevD+Befjnbd3d3BwA0adJEluIiMzMT2dnZeOKJJwz2EbmfE2dnZ5SVlVW5Xd3hzLp2JEkS5pQNJBaWwGQ1lixZApVKhSlTphj9GJVKBVdXV/1F903X2dkZrq6usp97RYC57AZsbGyEyywq3fMsWs+CJEk8RxEJRax3GFVbR48exddff4327dsrHaVcuj/sohUBNjY2wpzwT3S655mFC5F5ifUOo2opLy8Po0aNwurVq2WblCo3UQsXlUrFwsVCdM+zaEUACxcSDQsXUlxERAQGDRpk1YfU6r5Fi1a4cKjIckQeKhItM9VsnJxLitq0aRNOnDiBo0ePKh3lsUTtceFQkeVwqIjIMli4kGJu3LiBt99+G3v27NEf7mmtRC1cOFRkObp9Q7QigIULiYaFCynm+PHjuH37Np588kn9bWVlZThw4AC+/PJLFBUVwdbWVsGEfxO1cOFQkeWIXLiI1ktENRsLF1LMP/7xD5w5c8bgtvHjx6Nly5aYOXOm1RQtgLiFC3tcLEer1UKlUglZuIiWmWo2Fi6kGBcXl4fOIVS7dm00aNDA6s4tJHLhIlpmUekKF9GwcCHRsH+QyAi6P+yi9V5wqMhyRB1yEbXgopqLPS5kVaz1rLGi9rjwqCLL0Wq1QhYu7HEh0Yj3LiNSgKiFC+e4WI6oPRcsXEg0LFyIjMAF6KqHsrIys7Ut6lCRqLmp5uLeSmQEUXtcOFR0j1qtxty5c/Hiiy/iwoULZtmGOYaKcnNzMW3aNHzyyScAgH//+9+ytg+wx4XEwzkuREYQtXDhUBFw4cIFPPPMM8jMzIRWq0V6erpZtmOOoaLs7Gx8/fXX+v2uuLhY1vYBFi4kHva4EBlB1MKFQ0X3XrPS0lL9a6jRaMy2Hbl7XPz8/DBp0iT99W7dusnaPiDePk3EwoXICCIXLjW9x6VVq1Y4deoUhg0bBgA4e/asWbZjrqOKFi1aBGdnZwBAixYtZG+fPS4kGhYuREYwZ+GSkJAge5s65hoqkiQJ69atk71dc2nQoAG+/fZbDB48GFevXjXLNsx1VJGzs7O+p8Vck2hZuJBIWLgQGcGchcuePXtkb1PHHENFWq0WUVFR+PPPP2Vt1xLGjRuHsWPHmqVtcx6dM2DAALO0C7DHhcTDwoXICOYqXM6fP49z584BAAoKCmRtG5B/qOj8+fPo06cPvvnmG9SqVUu2di3Fzs7OLPNEAPMWAPXr1zdLuwALFxIPCxciI5ijcCktLcWrr76qb/P48eOyta0j51BRSUkJ+vXrp8/ZoEEDWdqtLrgAHZFlsHAhMoI5zlV08uRJxMfH6wuXEydOyNa2jpwfSPb29vjnP/+JevXqAQBcXV1la7s6EHUhN1ELLqq5xHuXESnAHH/Yu3Tpgv379yM8PBwA0KdPH9m3IfdQ0csvv4yEhARMmzYNPXv2lK3d6kDkcxWJmJtqLu6tREYw1xyXzp076+dctG/fXta2AfMcDl23bl289957CA4OlrVd0YlcuLDHhUQi3ruMSAGiruOiUqmEyywqUQsAUXNTzcXChcgIohYuXIDOctjjQmQZ4r3LiBQgcuEiWmZRiTpXhIULiUa8dxmRAkQtXHiSRcsR9egcFi4kmkoVLitWrECTJk3g5OSE4OBgHDly5LH3z87ORkREBLy9veHo6IiAgADs3LmzUoGJlCBq4cKhIsvhUBGRZdiZ+oDNmzcjKioKq1atQnBwMD799FOEhYXh4sWL8PDweOj+xcXF6N+/Pzw8PPCf//wHjRo1wvXr11G3bl058hNZhDnWcbEEDhVZjqhDRaL2FFHNZXLhsnz5coSHh2P8+PEAgFWrVmHHjh1Ys2YNZs2a9dD916xZg7t37+Lw4cOwt7cHADRp0qRqqYksTNQeFw4VWY6oBQB7XEg0Jn09KC4uxvHjxxEaGvp3AzY2CA0NRWxsbLmP+fXXXxESEoKIiAh4enqibdu2+PDDD1FWVvbI7RQVFSE3N9fgQqQkUQsXDhVZDoeKiCzDpHdZRkYGysrK4OnpaXC7p6cn0tLSyn3M1atX8Z///AdlZWXYuXMn5s2bh08++QQffPDBI7cTHR0NNzc3/cXX19eUmESyUqvV+P333wEAMTExuH37tsKJjKPValFUVAStVsvi38yKi4tRUFAArVaLwsJCpeMYRZIk7Nu3D1lZWbh27RpOnz6tdCQio5j964FWq4WHhwe++eYbdOrUCcOHD8fcuXOxatWqRz5m9uzZyMnJ0V9u3Lhh7phEj/TDDz9g+vTpAIAPP/wQK1asUDiRcd544w3MmDEDkiTBx8cHf/75p9KRKnTy5El4eHigpKQEUVFReOutt5SOZJRevXrhu+++Q2JiIho1avTIL3LWJD09Hc8++yxu3LiBvXv34uWXX1Y6EpFRTCpc3N3dYWtri/T0dIPb09PT4eXlVe5jvL29ERAQAFtbW/1trVq1QlpaGoqLi8t9jKOjI1xdXQ0uREp5/vnn9fOzAGDEiBEKpjFehw4d9D87ODigZcuWyoUxkpeXF0pKSgDc+9Lj7e2tcCLjdOrUSf+zl5cX3N3dFUxjHC8vL/Tt21d/fdy4ccqFITKBSYWLg4MDOnXqhH379ulv02q12LdvH0JCQsp9TI8ePXD58mWDcfbExER4e3vDwcGhkrGJLMfT0xMjR44EAPTu3VuIAgAAJkyYABcXFwDA66+/joYNGyqcqGLe3t549dVXAdz7AjNp0iSFExnnnXfe0c8TmT9/PuzsTD7uQRHz5s0DcO9vuyjPNZHJQ0VRUVFYvXo11q9fj/Pnz2PSpEnQaDT6o4zGjBmD2bNn6+8/adIk3L17F2+//TYSExOxY8cOfPjhh4iIiJDvtyAys9mzZ8PT0xPz589XOorRnJycMG7cOKhUKmGGXIB7f2NUKhUGDx4szLIJTzzxBDp16gQXFxe8+OKLSscxWnBwMDp16oTw8HB9kUtk7Uz+WjB8+HDcuXMH8+fPR1paGjp06IBdu3bpJ+wmJycbzKz39fXF7t27MXXqVLRv3x6NGjXC22+/jZkzZ8r3WxCZma+vL65evap0DJNFR0dj4cKFcHJyUjqK0by9vZGQkIDGjRsrHcUkf/zxB0pKSgyGxUVw4MABpSMQmaRS/ZmRkZGIjIws9/9iYmIeui0kJARxcXGV2RQRVYFKpRKqaNERca0nGxsbODo6Kh2DqNoTb9EBIiISXmZmJgAgISEBwL1lB4iMIcYMMiIT/f7777h58yaKi4uRmJio/yNpjXS9kXFxcbh79y7UajUuXryocKqK6XKfOHECmzZtUjiNcc6ePYusrCwAsPpe4Ly8PKSkpODGjRsG+4i1Sk1NBXDvkHadK1eu6H++du0aAODQoUMAgEuXLgEA6tSpg4yMDNSqVctSUUlwKkmApUBzc3Ph5uaGnJwcHhpNBiRJMvim9vPPP+snihORdVi4cCGefvpp3L59Gx06dEB8fDz69OmDmJgY/b+hoaGoXbu20lFJZub4/GaPCwlNrVbDzc3tkf8fHR2N27dvo3nz5lXajiRJyMjIMMshxUeOHMH69esxduxYBAYGQqPRwMfHR5a2CwoKUFpaapYjRnS5hw0bhqeeekrWtu/evQs3NzfZJ7pevnwZd+7cwQ8//ICxY8eia9eusrWt1WqRmZkp2z6Sn5+PO3fuwN/f32AfkTOzbjtarRZ16tSpUjt37twBAP3vL0kSrl+/jqVLl8LFxQWtW7fG7du3hTlSjKwXCxdS1MqVK7Fy5UokJSUBANq0aYP58+djwIABRj3excUFOTk5+uuHDh3CwIEDsXXrVnTr1g12dnY4dOgQ+vfvX6WcBw8exJQpU5CUlIR69epVqa0HOTk5Yf369ejduzd69OiBjIwMdOzYUZa2p0+fjkuXLuGXX36Rpb376XJ36dJF9sXLfH198fHHH8u+2N/evXuRkpKCH374Ab1795Z1tdhdu3Zh+vTpuHXrlizDHpmZmbhw4QJ69OhhsI/IvcLtpEmTkJubiw0bNlSpHd3wZmBgIIB7hcuxY8ewdOnSKmckuh8LF1JU48aNsWTJErRo0QKSJGH9+vV47rnncPLkSbRp06bCx6tUKoPuR11Xs7OzM1xdXVFQUCBLzry8PJSWlgpzHhodjUYj5HmK1Go18vLylI5hktzcXP25oUSSm5sr3HNNNRsLF1LU4MGDDa4vXrwYK1euRFxcnFGFi6XoVkUV7UOJZ4e2HN3zLNoZorVaLc8OTUJh4UJWo6ysDFu2bIFGo3nkKSSUovvDLsBcdgMsXCxH1MJFkiQWLiQUFi6kuDNnziAkJASFhYWoU6cOtm7ditatWysdy4Duw4iFCz2KyIWLaJmpZuPeSooLDAxEfHw8/vrrL0yaNAljx47FuXPnlI5lgENFVBFRCxcOFZFo2ONCinNwcNAfrtypUyccPXoUn332Gb7++muFk/2NPS5UEd2+IVrhwh4XEg33VrI6Wq0WRUVFSscwIGqPi0qlEi6zqHTPs2i9F+xxIdGwx4UUNXv2bAwYMAB+fn5Qq9XYuHEjYmJisHv3bqWjGRB5cq5omUWl1WqF7Lng5FwSDQsXUtTt27cxZswYpKamws3NDe3bt8fu3burvGCc3DhURBURuXARMTfVXCxcSFHfffed0hGMIupQEQsXyxG1cOFQEYlGvHcZkQLY40IVEbVwYY8LiYZ7K5ER2ONCFRG1cGGPC4lGvHcZkQJEnpzLwsUyRC1cODmXRCPeu4xIAaIOFdna2rJwsRCtVgtbW1ulY5iMhQuJhoULkRE4VEQVsUSPiznO4szChUTDwoXICKL2uLBwsRxLFC7Tp0/HwoULkZmZKVubnJxLouHeSkb5/vvv0aBBg4dWtB0yZAhGjx6tUCrLYY9L9SNJEn788UecP39elvYsMcm1pKQES5cuxciRI2Vrk5NzSTQsXMgoL730EsrKyvDrr7/qb7t9+zZ27NiBCRMmKJjMMjg5t/rZsGEDxo8fj0WLFsnSnqV6LurXr4+5c+fK1h6Hikg0LFzIKLVq1cLIkSOxdu1a/W0//PAD/Pz80KdPH+WCWYg5CxdzFhYqlcpsxVZ+fj6Ae/Mu/ve//2HZsmU4deqUWbYlp+TkZADAH3/8AQD4/fffUVpaWuV2zTlUlJWVBQDw8fHB0aNH0bt3b9naZuFCouHKuWS08PBwdOnSBbdu3UKjRo2wbt06jBs3rkb80TNn4RIXFwcAyM3Nlb1tGxsblJWVyd4uAH0Re3+PRV5eHoKCgsyyPTkUFRXho48+MrhNo9EgISEBHTp0qFLb5hxy0Q1nzZo1C15eXrK2Leph3FRzcW8lo3Xs2BFBQUH4/vvvcfz4cZw9exbjxo1TOpZFmHNybps2bQAAS5YsQWpqqqxtm/Mki2PHjgVgeDbknj17mmVbcnF0dNQXLkuXLsXmzZvx8ccfo0WLFlVu25yHQ4eEhAAAateuLXvbnJxLouHeSiZ59dVXsW7dOqxduxahoaHw9fVVOpJFmHNyrpubGwAgPT0dgwcPxpEjR2Rr25xzXOrUqQPg3j5Rq1Yt2Nvbo0ePHmbZlpx0uevXr4//+7//Q0REhCwFgTl7LszZq8nJuSQaFi5kkpEjR+LmzZtYvXp1jZiUq2Opw6Hz8vLw/fffy9aeJSbnBgQE4MCBA9i0aRNq1apl1m1ZM1GHXNjjQqLhHBcyiZubG1544QXs2LEDQ4YMUTqOxVjiqCInJye89tprGDhwoGxtWuqootatW6N169Zm3441E7lwYY8LiYSFC5ns1q1bGDVqFBwdHZWOYjG6DyRzFgFLly5F//79kZGRIVubPBzackQtXETNTTUXCxcyWlZWFmJiYhATE4OvvvpK6TgWZYkeF2dnZ9nbZOFiOaIWAOxxIdGwcCGjdezYEVlZWfjoo48QGBiodByL4sq5VBFRCxdOziXRsHAhoyUlJSkdQTE8VxFVRNQCgJNzSTTcW4mMwCX/qSKi9riItk8TifcuI1KAJSbnmoOoQ1wiErlwETE31VzcW4mMIHKPC8DCxRJELlxEHOKimku8dxmRAkQtAETNLSJRCxdRc1PNxb2VyAjscaGKiFoAsMeFRCPeu4xIAVlZWQCAjIwMs51t2RxYuFiOiIXL3bt3UVJSgvz8fBQWFiodh8goYr3LiBTw448/olu3bgDunVBw3rx5CicyTlFRkf5s0+fPnxeut0gk6enpyMzMRGFhoexn+DaXpKQkNGnSBElJSVi9ejX69++vdCQio7BwIapAp06dDL5Ji3AGZACIiorC3LlzAQC9e/fG1q1bFU5UsTNnzqBFixYoKSnB3LlzMWfOHKUjGaVbt27Yvn07zpw5g8DAQOTm5iodqUKNGzeGj4+P/roo+zURCxeiCjRr1gzPP/88AMDX11fWkyCaU1hYmP5nR0dH9OzZU8E0xqlbty7S09MBABqNBnZ2YqyR+cwzz+h/Dg4Ohqurq4JpjGNnZ4fZs2cDAGxtbTF16lSFExEZh4ULkRFmzZoFlUqFN998U5iJjIMHD4avry8AYMyYMfDw8FA4UcV8fX0xZswYAPc+WN966y2FExnnnXfe0f+8cOFC5YKYaOTIkXB2dkZYWBg8PT2VjkNkFBYuREZo1aoVcnJyEBERoXQUo6lUKsyYMQP29vaYMmWK0nGMNmPGDADA008/DXd3d4XTGKdp06bo1q0bmjVrJtSQi729PdLS0rBlyxaloxAZTYx+WCIrYGtrq3QEk02YMAETJkxQOoZJ/Pz8EBMTg/bt2ysdxST79u1TOkKliLhfU83GwoWEJkkS1Gq1/rru5zt37uDmzZuwt7cHAOTk5CiSzxj5+fn6fwsKClBSUmLVeXV0uYuLi2XPGxAQgMLCQtkP0ZUkSd9mfn6+VT/PGo0GZWVlyMnJMdhHrDVzUVERgL/fa5IkQaPRKBmJqikWLiQ0tVoNNze3h25/5ZVXAABvvfUWbG1tcfjwYUtHM9rZs2cBALt27UJSUhIkScJvv/2mcKqK6XIfOnTIaj9MH1RSUoKLFy8CuPd837hxQ+FEj6bVaqHVavHHH38Y7CPWmrm0tBQADCZU339oeEJCAgAgMzPTssGo2lFJAizukJubCzc3N+Tk5AgxW58s58Eel+3bt2PUqFEKJiKiBy1cuBBDhgxBUlISnJycUFhYCHt7e5SUlMDf3x/Xr1/HU089hbp16yodlWRmjs9v9riQoqKjo/Hzzz/jwoULqFWrFrp3746PPvoIgYGBRj1epVIZvBkaN24MAJg5cyYCAgJQXFyMxMREtG3b1iz55RAXF4fVq1cjPDwcAQEBUKvV8Pf3VzpWhXS5R4wYgdDQUKXjGOXs2bPIysrC2rVrER4erl9Y0Brl5eUhJSUFAQEBBvuItWbW9a54e3vrb7t69SoWL14MFxcX/WJ33bp1Q0xMDNq0aYP4+HgUFBQAAFfuJaOxcCFF7d+/HxEREejSpQtKS0sxZ84cPP300zh37hxq165d6Xb79euHnj17oqCgAAcPHqzyqqA5OTlYt24d3nrrLbMcDr169Wp069YNPXr0QEZGBjp27ChLu3FxccjNzcXTTz8tS3sPWr16NZ588km8/PLLsra7du1ahIWFGSyQJoe9e/ciJSUFa9euRbdu3WTNnZ2dje+//162Q+YzMzNx4cIF/VFKun1E7uc6NjYWGo2mysWnbghO96VDkiQcO3YMixcvLvf+ut6V1q1b4/bt23BxcanS9qnm4OHQpKhdu3Zh3LhxaNOmDYKCgrBu3TokJyfj+PHjSkczcPDgQcyZMwd37txROopJvvnmGyxfvlzpGCabOnUqduzYoXQMkxw4cACzZ8/Wn9dKFKtWrcI///lPpWMQGY2FC1kV3STP+vXrK5yEyDQCTBckqhZYuJDV0Gq1mDJlCnr06GHVc1KIiEg5nONCViMiIgIJCQk4ePCg0lEeopuzIOK3ahEzi0jU51mSJGFOY0EEsHAhKxEZGYnt27fjwIED+iODqOr4gWR5fM6JzIuFCylKkiS8+eab2Lp1K2JiYvDEE08oHYmIiKwYCxdSVEREBDZu3IhffvkFLi4uSEtLAwC4ubmhVq1aCqf7G4eKqCKiPs8cKiLRcHIuKWrlypXIyclBnz594O3trb9s3rxZ6WjVAj+QLI/POZF5sceFFCXqt1QiIlIGe1yIjCDyUBHR43CoiETDwoWommOxZRl8noksg4ULUTXGb9KWx+ecyLxYuBAR1WAcKiLRVKpwWbFiBZo0aQInJycEBwfjyJEjRj1u06ZNUKlUGDJkSGU2S0SVwCEMy+DzTGQZJhcumzdvRlRUFBYsWIATJ04gKCgIYWFhuH379mMfl5SUhOnTp6NXr16VDkukFFG/kYqaW2QiPuciZqaay+TCZfny5QgPD8f48ePRunVrrFq1Cs7OzlizZs0jH1NWVoZRo0Zh0aJFaNq0aYXbKCoqQm5ursGFiIiIyKTCpbi4GMePH0doaOjfDdjYIDQ0FLGxsY983HvvvQcPDw9MnDjRqO1ER0fDzc1Nf/H19TUlJpHZmHs4wBztcwjDMkR9nkXNTTWXSYVLRkYGysrK4OnpaXC7p6enfqn2Bx08eBDfffcdVq9ebfR2Zs+ejZycHP3lxo0bpsQkkp0lutIvXbqEXr164dSpU7K1aYncRUVFZt+GSEQcdhExM9VcZj2qSK1WY/To0Vi9ejXc3d2NfpyjoyNcXV0NLkTV3WeffYaMjAxs2rRJ6ShG0Wg0AIB58+bh8uXLCqcxXn5+PgBg/fr1KCgoUDgNEZnKpCX/3d3dYWtri/T0dIPb09PT4eXl9dD9r1y5gqSkJAwePFh/m1arvbdhOztcvHgRzZo1q0xuIkWYo1tdN8yq+xA9ePAgbt68icaNG8u+LbmsW7cO06dPB3DvOTl79iyaN2+ucKrHKykpweTJk5GSkgIAOHDgAG7evIkWLVoonKxi586dAwD897//RXx8PIqLizFz5syHer8rg0NFJBqTChcHBwd06tQJ+/bt0x/SrNVqsW/fPkRGRj50/5YtW+LMmTMGt7377rtQq9X47LPPOHeFhGHOrvSff/4ZAODn5wetVqufnC4Xc3wwrVq1CsXFxfrrSUlJsm9DbpIkQa1WGzwfcv4NMmcBcO3aNQDAtm3bYGNjA61Wi7Fjx8pSuAAcKiKxmDxUFBUVhdWrV2P9+vU4f/48Jk2aBI1Gg/HjxwMAxowZg9mzZwMAnJyc0LZtW4NL3bp14eLigrZt28LBwUHe34ZIQEuXLgUALFiwAHv37sXWrVvRunVrWdo21wdSbGws5s6dCwBwc3NDvXr1zLIdOTk4OGDdunV49tlnAQCurq5wcnKSfTvmeM4HDRqk/1mr1aJdu3Zo37697NshEoHJhcvw4cOxbNkyzJ8/Hx06dEB8fDx27dqlr/yTk5ORmpoqe1Aia8Bu9XtUKpV+aYN3330XY8aMUTiRcWxtbTFgwAAAQEREhMJpTPfSSy8BACZPnixbgcR9mkRj0lCRTmRkZLlDQwAQExPz2MeuW7euMpskUpTIXen8YCqf3HNyLPE8h4WFYcqUKQgKCpK1XZH3b6p5KlW4EJEY+IFkeeZ8zlUqFTp06GC29olEwJMsEpmAvRdU3XCfJtGwcCEyAnsuqDrj/k0iYeFCVI2pVCp+o7YQ3fPMIoDIvFi4EJmARQBVN9ynSTQsXIiMwG/RVJ1x/yaRsHAhqsb4gWR5fM6JzIuFC5EJROxWFzGziER9nkXNTTUXCxciI/BbNFVn3L9JJCxciKoxfiBZHp9zIvNi4UJUgc2bN2PSpEkAgGeffRb79+9XOBFR1UmSpN+fY2JiMGTIEKUjERmFhQtRBXJycnDjxg0AwOXLl1FUVKRwIuNcuXIFKSkpyMnJQXx8vNJxqrWMjAycO3cOAHD48GGUlZUpnMg4V65cgUajgUajwZUrVzjfhYTAwoWoAmPGjEHdunUBAC1btkT//v2VDWSksLAw7N69G5cuXUKPHj2QmZmpdKQKnTx5Ep6enigpKcG0adPw1ltvKR3JKDNmzMAnn3wCABg8eDB27NihcKKKqVQqzJ07V3/93Xff5TAXCYGFC1EFnJycMGbMGADA1KlThfnjPnz4cP3PTz31FBo0aKBgGuN4eXnpe7TKysrg5eWlcCLjvPTSS/qfXV1d0bdvXwXTGG/YsGFwdnZG7dq18eKLLyodh8goLFyIjPDBBx9g1apVGDVqlNJRjDZlyhTY2toCAObNm6dwGuN4e3vj1VdfBQA4ODhg8uTJCicyzjPPPANfX18A94pbFxcXhRMZx87ODrt27cKuXbv0+wqRtbNTOgBRVUiSBLVarb+em5sL4N7YvYuLCxwcHPTXq6p79+64evVqldt50O3bt/X/3rlzB/n5+bLkBYCnn34aly5dgoeHh2xt6uhy5+bmytr2Cy+8gJUrV6Jv377IzMyUdYhLkiRkZGQAuJdfztwTJkzAP//5T4SFhcnSbn5+PoqLi3HlyhWDfUTu11E3DFrVdrOzsw3akSQJqampVWqTqDwsXEhoarUabm5uD92u+6a+d+9eAEBhYaFFc5mipKRE/29paSm0Wq1seRctWgTAPL+/LreceYF7Qy0//fQTvL29zZJbN3G2pKRE1vYHDRqEQYMGAZDn+S4pKYEkSSgsLDTYR6x1X9Y9r/fny8vLe+j/iaqKhQsp6sCBA1i6dCmOHz+O1NRUbN261aTDMl1cXJCTk6O/fujQIQwcOBBbt25Ft27dYGdnh0OHDqFNmzZmSC+PM2fOAAAaNWoEb29vZGRkWHVeHV3uunXryp7XXL9/amoqPD09Adx7vq35ec7MzERBQQHatGljsI9Ya+aLFy8CAAIDAwHc63HJz8/X//+RI0f0txNVBQsXUpRGo0FQUBAmTJiA559/3uTHq1QquLq66q/Xrl0bABAXF4e7d++iuLgYiYmJVn1ETVxcnP7fu3fvQq1W6z8ErJku94kTJ7Bp0yaF0xjn7NmzyMrKAvB3fmuVl5eHlJQU3Lhxw2AfsVa6YaGTJ0/qb7t/aNXPzw+XLl3CqVOnLJ6NqheVJED5m5ubCzc3N+Tk5Bh8SFH1olKpTO5xedAvv/zChbSIrMzChQsxZcoUxMTEoHHjxrh58ya8vb2RmpqKPn36ICYmBqGhofovHlR9mOPzmz0uVK3Uq1cPADBz5kwEBAToe1zatm2rcLJHi4uLw+rVqxEeHo6AgACo1Wr4+/srHatCutwjRoxAaGio0nGMoutxWbt2LcLDw9GtWzelIz2SrsclICDAYB+x1sy6Hhdvb2/9bVevXsXixYsNjrJq3rw5bt68qZ8Lk5ycbNmgJDwWLlQt9evXDz179kRBQQEOHjxY5UXjjh07hnnz5mHnzp1mWcdl9erV6NatG3r06IGMjAx07NhRlnbfe+89uLi4YOrUqbK096DVq1fjySefxMsvvyxru88++yxmzJiBnj17ytru3r17kZKSgrVr16Jbt26y5k5NTcULL7yAnTt36o/UqYrMzExcuHABPXr0APD3PiL3c71+/XqcOnUKy5cvr1I75c1xOXbsGBYvXlzu/YOCghATE6MveNRqNXtcyChcx4XICAkJCThw4IDSMUwWGxurn9gpkgMHDuD8+fNKxzDJrVu3cOrUKdy8eVPpKCY5c+YMDh48qNj2u3TpAgBwdHRULAOJhYULkRGKi4vh4OAgzKq5OoWFhfxAsBDdmkGinMtKx8HBQX+4tRJsbGz0OYiMwcKFyAi6wkU0RUVFcHJyUjpGjaB7nkUsXIqLi5WOQWQ0znEhReXl5eHy5cv669euXUN8fDzq168PPz8/BZMZErlwETG3iHQ9W9a6QNyjsHAh0bBwIUUdO3bM4IR0UVFRAICxY8di3bp1CqV6WElJCezsxHu7sMfFcnSFi2g9Lvb29ixcSCji/SWmaqVPnz5CrKQpao8L57hYjsiFi5JzXIhMxTkuREYQtXApLi5m4WIhup4tDhURmRcLFyIjiFq4sMfFckTtcWHhQqJh4UJkhNLSUtjb2ysdw2Sc42I5NjY2sLOzE65wsbe3R1lZGbRardJRiIzCwoXICCL2uJSVlaG0tFS43CJzcnISrnDR7R+c50KiYOFCZAQRCxfdByh7XCzH0dFRyDkuADhcRMJg4UJkhOLiYuGGinQfoJzjYjmOjo7CFQAsXEg0LFyIjFBSUiJsjwsLF8sRscdFV5CzcCFRsHAhMoLIQ0UsXCyHc1yIzI+FC5ERRC5cOMfFchwcHIQtXNjjQqJg4UJkBM5xIWOI3OPCwoVEwcKFyAgiznHRfRCxcLEcznEhMj8WLkRGMOdQkVqtNku75uxxEeH8Ug9KTU3Fn3/+iYyMDLNtQ+6jin777TfUrl0biYmJsrX5IM5xIdGwcCEygrmGii5duqQ/I7bcxYC55rj8+uuv8Pf3x7lz52Rt11ymTJkCHx8fNG/eHMuWLcPBgwfNti25e1x0RcVHH30EANiwYYPsQ1HscSHRsHAhMkJJSQns7OQ/mfqmTZv0S63//PPPsrat+4CTq6eotLQUERERGDFiBDIzM3H69GlZ2jW35ORk5OTk6K/7+PiYbVtOTk6yFi7/+Mc/0KNHD/312NhY2QsX3f4h2twcqrlYuBAZwRxnWdZqtVi7dq3++m+//YbMzEzZ2pe7x+X999/HunXr9NcvX74sS7vm9vPPP+OLL76Ara0tAKBJkyZm25Y5Tlj4yiuv6H+OioqCq6urrO1zqIhEw8KFyAjmmOMyf/58pKeno0GDBgCAOXPm6H+Wg9xzXN555x3MmTMHw4YNg4eHh1AfdBMmTMCePXswYsQINGzY0GzbMcdRRUOHDkWHDh0AAE2bNpW1beDvoSKRXk+q2eTv+yaqhsxRuIwYMQL+/v5wcXHBxIkTZf9Q0n3zlyt3nTp1MHfuXAD35uP8+9//Rnh4uCxtW0JwcDDUajVSUlLMtg1zHFXk4uKCN998ExMnTpS1XR0eDk2iYY8LkRFKS0tln5zbpk0bs37wFxYWwtHRESqVSva2VSoVbGz45+NBIp+riD0uJAr+5SEygqgr53LVXMsScR0X9riQaFi4EBlB1MJFtMyiE3HlXFtbW9jY2LDHhYTBwoWoApIkCVu4sMfFskQ8VxFgnqOhiMyFhQtRBcrKyiBJkpDnKuJy/5YlYo8LwMKFxMLChagCch+dYykcKrI8R0dHFBUVCXdKBBYuJBIWLkQV0P1BF63HhUNFlqfr4RKtCLCzs+McFxIGCxeiCojc48KhIsvSPd+iDRexx4VEwsKFqAK6b6Ii9riwcLEsXQ+XiIdEs3AhUbBwIaoAe1zIWKKuiWJvby9cZqq5WLgQVUDUwqWwsJBzXCxM5B4XznEhUbBwIXqM/fv3Y/LkyQCADz74AHFxcQonMh6PKrI80ea4FBcXY+rUqbh+/Tp27tyJd955R+lIRBVi4UL0GElJSfpiZf/+/Th37pzCiSqWn5+PoUOH4sSJE/jzzz8xb948pSMZ7dSpU/Dz80NJSQlmzpyJ6dOnKx3JaLGxsZg9ezYAYNy4cdi5c6fCiSqm1WqxceNGZGdn49atW/j999+VjkRUIRYuRI8xfPhwuLq6AgAaNmyIUaNGKZyoYiqVCrGxscjLy0NaWhr27t2rdCSjubu7IycnB8C9Xgvdcy+C+4vcc+fOIT09XeFEFXNycjIoDqdNm6ZgGiLj2CkdgKgqJEmCWq3WX8/KygIAnDhxAsXFxfo5B3/99Velt/HUU0/hv//9L4YOHYr4+Pgq5S3PlStX9P82atQIpaWlVcoLAM899xx++OEHAMDIkSOr3F55dLnv3Lkja/sDBw7Er7/+Cjs7O4SEhMjatiRJSE5OBnAvv5xt+/v7w9PTE+np6XB2dkazZs2q1H5paSkKCwvx119/Gewjcr+WXbp0gZ2dHWxsbODv71/p9nXzerKzs/W3Xbp0SY6IRAZYuJDQ1Go13NzcHrpd12V/7NgxaDQaNGnSpNLb+PDDD+Hh4YHp06fDzk7+t0zDhg31/zZo0AC5ublVygsAb731Fn744Qd4e3vjpZdegkqlkiGpIV1uFxeXKue934wZM/Drr7+iX79+aN++vWztAsDp06dRv359APfyy5kbuNdjMWPGDAwcOBCBgYFVakutVuPWrVto0qSJwT4id2YA+Oc//wlbW1s0b9680m2kpKQAAHx8fB66jUhOLFxIcStWrMDSpUuRlpaGoKAgfPHFF+jatatRj3VxcdEPLQDAoUOHMHDgQGzduhXdunWDra0tEhMT4eHhUaUP788//7zSj62IbjjE1dUVderUwd27d6uc19PTE6tWrULHjh3h5eUlV1QD9w/jVDXv/Tw9PREXF4eWLVuaZe0cXfHp6uoKT09PWduePHkyNBoNJk+ejDp16lSpLVtbW1y9ehUeHh4G+4jcmQFgwoQJVW4jKysLGo1Gvy9IkoTGjRvr///+nhiiquAcF1LU5s2bERUVhQULFuDEiRMICgpCWFgYbt++bdTjVSoVXF1d9ZfatWsDAJydneHq6gqNRoORI0caDCdZM41Gg7CwMFnyjh49Gm3btpUh1ePNmzdP9ue3Xbt2ZilaCgoKMHXqVNnb1VGpVJgxY0aVixYAyMvLw9ChQ4Xad3v37v3IvKdOnbJwIqquKtXjYso35NWrV+P7779HQkICAKBTp0748MMPjf5GTdXb8uXLER4ejvHjxwMAVq1ahR07dmDNmjWYNWuWye3pTm6n0WiQk5OD3NxcAEBubq7VnvguPz9f/29eXh4A686ro8sNiJEXMMycn59v0FtnbXQFQG5ursE+Yq2ZNRoNgL/3BUmS9PuzJEkICAjAxYsXkZmZqWRMqg4kE23atElycHCQ1qxZI509e1YKDw+X6tatK6Wnp5d7/5EjR0orVqyQTp48KZ0/f14aN26c5ObmJt28edPobebk5EgApJycHFPjkhUrKiqSbG1tpa1btxrcPmbMGOnZZ5+tVJs7d+6UAPDCCy9WdFm8eLGUnZ0tbdu2Tfrll1+kbdu26a/r/s3Ly5PhrwpZG3N8fps8VHT/N+TWrVtj1apVcHZ2xpo1a8q9/4YNGzB58mR06NABLVu2xLfffgutVot9+/Y9chtFRUXIzc01uFD1k5GRgbKysofG7D09PZGWllapNps2bSpHNCKSUfv27WFjc+/jpk2bNgBQ6fc4kUlDRcXFxTh+/Lj+iA0AsLGxQWhoKGJjY41qIz8/HyUlJfqZ/eWJjo7GokWLTIlGBAAIDAxEYmKivps9Ly8PTz31FPbv31+leQc5OTkoKSmBu7u7XFENZGZmokGDBrLl1UlNTUW9evXMtvT/jRs3MGTIENny6qSkpKBhw4ayz3PRPb/btm2Dr6+vrG0DQHJyMho1agRbW9sqt/XgvqDbR+RWWlqKW7duwd/fv0rtPGrfdXFxQbNmzVBQUAAA+vfQxYsXq7Q9qrlMKlwe9w35woULRrUxc+ZM+Pj4IDQ09JH3mT17NqKiovTXc3NzzfJHhpTl7u4OW1vbhxbqSk9Pr9KRMC1atND/rOut69ChQ5UWMxs/fjyys7OxdevWSrdhDLny6nh5eeH9999HeHh4ldsqj+7wWbny6vTv3x+fffYZXnnlFdnaBP5+fvv27Sv74nZarRY9e/bEunXrMGzYsCq3J/e+8Ch//PEHhg0bhrS0tCoVRhXldXBwQLNmzfQ9L0FBQYiPj8eNGzcqvU2qmSx6VNGSJUuwadMmbN269bHfAB0dHQ2OFBFp9UwynoODAzp16mQwbKgbRgwJCVEw2cMSExOF2w/z8vL0i6GJRq1WC3M0jY6NjQ1q1aqlXyxOFK6urigtLUVSUpJZt2Nvb4/WrVvrC5e6desCAK5du2bW7VL1Y1LhUpVvyMuWLcOSJUvwv//9T/ZFpUhcUVFRWL16NdavX4/z589j0qRJ0Gg0+qOMqkq3zouLi0uV2klOToafn58smR5HrrzAvWEiwHBBMLnJmdcSzJ3Xx8dHtkXXLPXc6vZr3YrClWVMXhsbm4d6Xjp27AgAQpwigayDSYVLZb8hf/zxx3j//fexa9cudO7cufJpqdoZPnw4li1bhvnz56NDhw6Ij4/Hrl27ZFtkS7fOS1UWRyspKUFKSopFChc58uroPkAbNWpU5bYeRc68lmDuvI0aNZKtcLHUc+vu7g4nJ6cqFy7G5n2w50VX6OhWByaqiMnruERFRWHs2LHo3Lkzunbtik8//dTgG/KYMWPQqFEjREdHAwA++ugjzJ8/Hxs3bkSTJk30M8nr1Kkj62Q+EldkZCQiIyOVjvFIt27dglartUjhIqdbt24BALy9vRVOUnP4+PgIcQbx+6lUKvj5+eH69esW22Z5PS+6f4kqYnLhMnz4cNy5cwfz589HWloaOnToYPANOTk52WAHXLlyJYqLi/Hiiy8atLNgwQIsXLiwaumJLED3TVS0wiUlJQUuLi7CDONUBz4+PkKdjVvHz8+vyj0uptL1vOiONiIyVqVWzn3cN+SYmBiD6+ae8EVkbiIXLuac30IP8/HxQWpqKrRarVA9CP7+/jh9+rTFt2tjY2OWE5dS9SbOO4tIIcnJyWjQoIH+PEiiYOFieT4+PigtLUVGRobSUUyiRI+Ljp2dHZo3by5UoUfK4p5CVAFLHVEkt5SUFLNOzKWH6Z5vuSboWoqfnx/S09NRWFho8W3b2tqiTZs2qFWrlsW3TWJi4UJUAVELl1u3brHHxcJ0z7duYrQodPv3zZs3FU5CVDEWLkQVuH79unCFiyRJHCpSgKenJ1QqlZA9LgAsemQRUWWxcCF6DEmShOxxyc7ORmFhIQsXC7Ozs4Onp6dwhUvjxo0BVH0ROiJLYOFC9BjZ2dnIy8sTrnDRfXCycLE8OVfPtRQnJyd4enqycCEhsHAhegyRD4UGzLtqLpVPztVzLUnJI4uITMHChegxRC1cuGqucnx8fISbnAuwcCFxsHAheozk5GTY29tXeBJRa5OSkoIGDRrA0dFR6Sg1johDRcC9RehYuJAIWLgQPcb169fh6+sr3OJYPKJIOT4+Prh9+zZKSkqUjmISXY+LJElKRyF6LLH+GhNVwooVK9CkSRM4OTkhODgYR44cMfqxljyi6MCBAxg8eDB8fHygUqmwbdu2Srdl7sIlOjoaXbp0gYuLCzw8PDBkyBBcvHjRbNuripUrV6J9+/ZwdXWFq6srQkJC8Ntvv5ltez4+PpAkCenp6VVua8mSJVCpVJgyZUrVg1XAz88PhYWFuHPnjtGPWbhwIVQqlcGlZcuWZkxJxMKFqrnNmzcjKioKCxYswIkTJxAUFISwsDDcvn3bqMdbsnDRaDQICgrCihUrqtyWuQuX/fv3IyIiAnFxcdizZw9KSkrw9NNPQ6PRmG2bldW4cWMsWbIEx48fx7Fjx9CvXz8899xzOHv2rFm2p5sQXdV5LkePHsXXX3+N9u3byxGrQrr93NThojZt2iA1NVV/OXjwoDniEemxcKFqbfny5QgPD8f48ePRunVrrFq1Cs7OzlizZo1Rj7dk4TJgwAB88MEHGDp0aJXbunXrllmPKNq1axfGjRuHNm3aICgoCOvWrUNycjKOHz9utm1W1uDBgzFw4EC0aNECAQEBWLx4MerUqYO4uDizbE9XMFZlnkteXh5GjRqF1atXo169enJFe6zKFi52dnbw8vLSX9zd3c0Rj0iPhQtVW8XFxTh+/DhCQ0P1t9nY2CA0NBSxsbEVPr6kpAQpKSnCHVGk1WqRmppq0TkuOTk5AID69etbbJuVUVZWhk2bNkGj0SAkJMQs22jQoAHs7e2rVLhERERg0KBBBvuuubm7u8PJycnkwuXSpUvw8fFB06ZNMWrUKE7wJbPj+cSp2srIyEBZWRk8PT0Nbvf09MSFCxcqfPytW7cgSZJwhcudO3dQVlZmscJFq9ViypQp6NGjB9q2bWuRbZrqzJkzCAkJQWFhIerUqYOtW7eidevWZtmWjY0NvL29K124bNq0CSdOnMDRo0dlTvZ4KpXK5EOig4ODsW7dOgQGBiI1NRWLFi1Cr169kJCQABcXFzOmpZqMhQvRI4i6houlV82NiIhAQkKCVc9tCAwMRHx8PHJycvCf//wHY8eOxf79+81WvFT2kOgbN27g7bffxp49e+Dk5GSGZI/n5+dn0vmKBgwYoP+5ffv2CA4Ohr+/P3788UdMnDjRHBGJWLhQ9eXu7g5bW9uHju5IT083al0W3R9w0QoX3aRQS6yaGxkZie3bt+PAgQP6891YIwcHBzRv3hwA0KlTJxw9ehSfffYZvv76a7Nsr1GjRpWanHv8+HHcvn0bTz75pP62srIyHDhwAF9++SWKiopga2srZ1QD/v7+OHXqVKUfX7duXQQEBODy5csypiIyxDkuVG05ODigU6dO2Ldvn/42rVaLffv2GTW/ITk5GQ0aNEDt2rXNGVN2KSkpsLGxgYeHh9m2IUkSIiMjsXXrVvz+++944oknzLYtc9BqtSgqKjJb+5XtcfnHP/6BM2fOID4+Xn/p3LkzRo0ahfj4eLMWLUDVV8/Ny8vDlStXuGIzmRV7XKhai4qKwtixY9G5c2d07doVn376KTQaDcaPH1/hYy19Vui8vDyDb6rXrl1DfHw86tevb1KOlJQUeHp6ws7OfG/viIgIbNy4Eb/88gtcXFyQlpYGAHBzc0OtWrXMtt3KmD17NgYMGAA/Pz+o1Wps3LgRMTEx2L17t9m2WdnCxcXF5aF5QrVr10aDBg0sMn/Iz88Pt2/fRkFBgVGv4/Tp0zF48GD4+/sjJSUFCxYsgK2tLUaMGGH2rFRzsXCham348OG4c+cO5s+fj7S0NHTo0AG7du16aMJueSxduBw7dgx9+/bVX4+KigIAjB07FuvWrTO6HUusmrty5UoAQJ8+fQxuX7t2LcaNG2fWbZvq9u3bGDNmDFJTU+Hm5ob27dtj9+7d6N+/v9m26ePjg6ysLKMLAGuh299v3ryJFi1aVHj/mzdvYsSIEcjMzETDhg3Rs2dPxMXFoWHDhuaOSjUYCxeq9iIjIxEZGWny45KTk/GPf/zDDInK16dPH1mWW7dE4SLSsvDfffedxbepe/5TU1PRtGnTKrUVExMjQyLj3L+WizGFy6ZNm8wdieghnONCVA5Jkize4yIXcy8+RxWTa/VcS9NNsOZaLGTNWLgQlSM7Oxt5eXnw9/dXOorJeIJF5cmxeq4SnJyc4OXlxcKFrBoLF6JyiLqGS0lJCW7fvs3CRWGurq5wdnYWrnABqn5kEZG5sXAhKoeoa7joju5h4aIslUpV6SOLlGbqInRElsbChagcycnJsLe3N+roI2uh1Wr1h1OLOsclJSUFr7/+OkpLS/H9999j8+bNSkeqtEaNGuHmzZvIz89XOopJ2ONC1o6FC1E5kpOT4evrCxsbcd4iERER6NevHwBg6NChOHHihMKJTHPnzh20a9cO33zzDSRJQlxcHNauXQsHBweDRQRFMGzYMMTFxWHTpk2oV68ecnNzlY5kNF3hItKRY1SziPNXmciCRDyiqFWrVvqfk5KS4OjoqGAa0zVs2BD/+te/oFKp9LedPXsWkZGRFj0sXQ5ZWVn6lXl9fHyEOuGgn58fioqKcOfOHaWjEJWLhQtROZKTk4U7omjkyJH6D/3p06ejTZs2Cicy3cCBA/Hcc88BuHfem3r16iE6OlrhVKZbuXKlfnn+iRMnGhRj1k6333O4iKwVCxeicojY4+Lu7o7AwEDUqlULCxYsUDpOpW3cuBGurq5Qq9XYsGGDcD1HANC8eXOMHTsWADBq1CiF05jm/kXoiKwRCxeiB5SUlCAlJUW4wgUA4uLicPXqVTg7OysdpdKuXLmC4uJiAPeGvET17bff4tKlS8KdgLJBgwaoVasWCxeyWlzyn+gBt27dgiRJQhYubm5ucHNzUzpGpRUXF+OVV17B8OHDERgYiFdffRVnzpwx65muzUWlUqF58+ZKxzCZSqXikUVk1djjQnSfP//8E4sWLQJwb00UrVarcKKaZe7cucjJycHnn3+OmTNnIiAgABMmTFA6Vo1y7do1ODo6Yu/evfj444/5HiCro5IEOOYtNzcXbm5uyMnJgaurq9JxqBp7/fXX8c033+iv//TTT3j++ecVTFRzxMTEoH///vjjjz/Qs2dPAPeGioKCgrBkyRJMmjRJ4YQ1Q6NGjfQL59na2iI7Oxt16tRROBWJyhyf3+xxIbrPiy++qP+5SZMmeOaZZxRMU7P06dMHJSUl+qIFuPca5OTksGixoLfeekv/c1hYGIsWsjosXIju069fPzg5OQEA1q9fL/QkV6LKmDZtmn5+18svv6xwGqKHsXAhuo+trS1ee+01DB8+HL1791Y6DpHF2dnZ4fvvv0dgYCCHSckqcY4LERERmQXnuBAREVGNxsKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiIRRqcJlxYoVaNKkCZycnBAcHIwjR4489v5btmxBy5Yt4eTkhHbt2mHnzp2VCktEREQ1m8mFy+bNmxEVFYUFCxbgxIkTCAoKQlhYGG7fvl3u/Q8fPowRI0Zg4sSJOHnyJIYMGYIhQ4YgISGhyuGJiIioZlFJkiSZ8oDg4GB06dIFX375JQBAq9XC19cXb775JmbNmvXQ/YcPHw6NRoPt27frb+vWrRs6dOiAVatWGbXN3NxcuLm5IScnB66urqbEJSIiIoWY4/PbzpQ7FxcX4/jx45g9e7b+NhsbG4SGhiI2Nrbcx8TGxiIqKsrgtrCwMGzbtu2R2ykqKkJRUZH+ek5ODoB7T0BNxGKNiEhsNfXzS/d7m9hH8lgmFS4ZGRkoKyuDp6enwe2enp64cOFCuY9JS0sr9/5paWmP3E50dDQWLVr00O2+vr6mxCUiIiIrkJmZCTc3N1naMqlwsZTZs2cb9NJkZ2fD398fycnJsv3iVDm5ubnw9fXFjRs32BOkML4W1oOvhXXh62E9cnJy4Ofnh/r168vWpkmFi7u7O2xtbZGenm5we3p6Ory8vMp9jJeXl0n3BwBHR0c4Ojo+dLubmxt3Qivh6urK18JK8LWwHnwtrAtfD+thYyPf6ismteTg4IBOnTph3759+tu0Wi327duHkJCQch8TEhJicH8A2LNnzyPvT0RERPQoJg8VRUVFYezYsejcuTO6du2KTz/9FBqNBuPHjwcAjBkzBo0aNUJ0dDQA4O2338ZTTz2FTz75BIMGDcKmTZtw7NgxfPPNN/L+JkRERFTtmVy4DB8+HHfu3MH8+fORlpaGDh06YNeuXfoJuMnJyQZdQt27d8fGjRvx7rvvYs6cOWjRogW2bduGtm3bGr1NR0dHLFiwoNzhI7IsvhbWg6+F9eBrYV34elgPc7wWJq/jQkRERKQUnquIiIiIhMHChYiIiITBwoWIiIiEwcKFiIiIhMHChYiIiIRhNYXLihUr0KRJEzg5OSE4OBhHjhx57P23bNmCli1bwsnJCe3atcPOnTstlLT6M+W1WL16NXr16oV69eqhXr16CA0NrfC1I+OZ+r7Q2bRpE1QqFYYMGWLegDWIqa9FdnY2IiIi4O3tDUdHRwQEBPDvlExMfS0+/fRTBAYGolatWvD19cXUqVNRWFhoobTV14EDBzB48GD4+PhApVI99uTJOjExMXjyySfh6OiI5s2bY926daZvWLICmzZtkhwcHKQ1a9ZIZ8+elcLDw6W6detK6enp5d7/0KFDkq2trfTxxx9L586dk959913J3t5eOnPmjIWTVz+mvhYjR46UVqxYIZ08eVI6f/68NG7cOMnNzU26efOmhZNXP6a+FjrXrl2TGjVqJPXq1Ut67rnnLBO2mjP1tSgqKpI6d+4sDRw4UDp48KB07do1KSYmRoqPj7dw8urH1Ndiw4YNkqOjo7Rhwwbp2rVr0u7duyVvb29p6tSpFk5e/ezcuVOaO3eu9PPPP0sApK1btz72/levXpWcnZ2lqKgo6dy5c9IXX3wh2draSrt27TJpu1ZRuHTt2lWKiIjQXy8rK5N8fHyk6Ojocu8/bNgwadCgQQa3BQcHS6+//rpZc9YEpr4WDyotLZVcXFyk9evXmytijVGZ16K0tFTq3r279O2330pjx45l4SITU1+LlStXSk2bNpWKi4stFbHGMPW1iIiIkPr162dwW1RUlNSjRw+z5qxpjClcZsyYIbVp08bgtuHDh0thYWEmbUvxoaLi4mIcP34coaGh+ttsbGwQGhqK2NjYch8TGxtrcH8ACAsLe+T9yTiVeS0elJ+fj5KSElnPBFoTVfa1eO+99+Dh4YGJEydaImaNUJnX4tdff0VISAgiIiLg6emJtm3b4sMPP0RZWZmlYldLlXktunfvjuPHj+uHk65evYqdO3di4MCBFslMf5Prs9vkJf/llpGRgbKyMv0pA3Q8PT1x4cKFch+TlpZW7v3T0tLMlrMmqMxr8aCZM2fCx8fnoZ2TTFOZ1+LgwYP47rvvEB8fb4GENUdlXourV6/i999/x6hRo7Bz505cvnwZkydPRklJCRYsWGCJ2NVSZV6LkSNHIiMjAz179oQkSSgtLcUbb7yBOXPmWCIy3edRn925ubkoKChArVq1jGpH8R4Xqj6WLFmCTZs2YevWrXByclI6To2iVqsxevRorF69Gu7u7krHqfG0Wi08PDzwzTffoFOnThg+fDjmzp2LVatWKR2txomJicGHH36Ir776CidOnMDPP/+MHTt24P3331c6GlWS4j0u7u7usLW1RXp6usHt6enp8PLyKvcxXl5eJt2fjFOZ10Jn2bJlWLJkCfbu3Yv27dubM2aNYOprceXKFSQlJWHw4MH627RaLQDAzs4OFy9eRLNmzcwbupqqzPvC29sb9vb2sLW11d/WqlUrpKWlobi4GA4ODmbNXF1V5rWYN28eRo8ejVdffRUA0K5dO2g0Grz22muYO3euwUmBybwe9dnt6upqdG8LYAU9Lg4ODujUqRP27dunv02r1WLfvn0ICQkp9zEhISEG9weAPXv2PPL+ZJzKvBYA8PHHH+P999/Hrl270LlzZ0tErfZMfS1atmyJM2fOID4+Xn959tln0bdvX8THx8PX19eS8auVyrwvevTogcuXL+uLRwBITEyEt7c3i5YqqMxrkZ+f/1BxoisoJZ5j2KJk++w2bd6weWzatElydHSU1q1bJ507d0567bXXpLp160ppaWmSJEnS6NGjpVmzZunvf+jQIcnOzk5atmyZdP78eWnBggU8HFompr4WS5YskRwcHKT//Oc/Umpqqv6iVquV+hWqDVNfiwfxqCL5mPpaJCcnSy4uLlJkZKR08eJFafv27ZKHh4f0wQcfKPUrVBumvhYLFiyQXFxcpH//+9/S1atXpf/9739Ss2bNpGHDhin1K1QbarVaOnnypHTy5EkJgLR8+XLp5MmT0vXr1yVJkqRZs2ZJo0eP1t9fdzj0O++8I50/f15asWKFuIdDS5IkffHFF5Kfn5/k4OAgde3aVYqLi9P/31NPPSWNHTvW4P4//vijFBAQIDk4OEht2rSRduzYYeHE1Zcpr4W/v78E4KHLggULLB+8GjL1fXE/Fi7yMvW1OHz4sBQcHCw5OjpKTZs2lRYvXiyVlpZaOHX1ZMprUVJSIi1cuFBq1qyZ5OTkJPn6+kqTJ0+WsrKyLB+8mvnjjz/K/fuve/7Hjh0rPfXUUw89pkOHDpKDg4PUtGlTae3atSZvVyVJ7CsjIiIiMSg+x4WIiIjIWCxciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBgsXIiIiEgYLFyIiIhIGCxciIiISBj/D/VtQEwIyFwmAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Stress discretization\n",
    "stress = data[pp.DISCRETIZATION_MATRICES][parameter_keyword][\n",
    "    mpsa_class.stress_matrix_key\n",
    "]\n",
    "# Discrete boundary conditions\n",
    "bound_stress = data[pp.DISCRETIZATION_MATRICES][parameter_keyword][\n",
    "    mpsa_class.bound_stress_matrix_key\n",
    "]\n",
    "\n",
    "\n",
    "T = stress * u + bound_stress * u_b\n",
    "\n",
    "T2d = np.reshape(T, (g.dim, -1), order=\"F\")\n",
    "u_b2d = np.reshape(u_b, (g.dim, -1), order=\"F\")\n",
    "assert np.allclose(np.abs(u_b2d[bound.is_neu]), np.abs(T2d[bound.is_neu]))\n",
    "\n",
    "T = np.vstack((T2d, np.zeros(g.num_faces)))\n",
    "pp.plot_grid(g, vector_value=T, figsize=(15, 12), alpha=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the traction on face i: T[2*i:2*i+g.dim] is the traction on the face as defined by the normal vectors `g.face_normals`. This means that for the bottom boundary, the traction T[bot] is the force from the interior of the box to the outside (since the normal vectors here are [0,1]), while on the top boundary, the traction T[top] is the force applied to to top faces from the outside (since the normals here point out of the domain).\n",
    "\n",
    "# What we have explored\n",
    "We have seen how linear elastic problems can be solved using the Mpsa class, which implements a multi-point finite volume method for momentum balance. The primary unknown in Mpsa is the cell center displacements, whereas tractions on faces can easily be post-processed."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
